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Abstract 


The primary focus of this study is the physics of the transition 
and early turbulence regimes in the time-developing mixing layer . 
Three-dimensional, time-dependent numerical simulations were carried 
out. In particular, we deal with the sensitivity of the mixing layer to 
the disturbance field of the initial condition. The growth of the 
momentum thickness , the mean velocity profile , the turbulence kinetic 
energy, the Reynolds stresses, the anisotropy tensor, and particle track 
pictures of computations are all examined in an effort to better under- 
stand the physics of these regimes. The amplitude, spectrum shape, and 
random phases of the initial disturbance field were varied. 

In carrying out this study, a new scheme of generating discrete 
orthogonal function expansions on some nonuniform grids was developed. 
In the present work it allowed us to compute in an infinite domain, 
eliminating image-flow problems . The new scheme retains the efficiency 
of the fast Fourier transform, but allows the application of more gen- 
eral boundary conditions by using a restricted set of mapping functions. 
In evaluations using linear test equations, the new scheme employed in 
the present work had errors up to six orders of magnitude smaller than 
in standard finite-difference methods (using equal numbers of grid 
points) . 

Due to computational limitations , all cases address the early or 
near field of the mixing layer. The results showed that large eddies 
may vary considerably, particularly in the turbulent structure measured 
by the anisotropy tensor. An interesting oscillatory behavior of the 
width of turbulent kinetic energy profile was observed . The most sig- 
nificant result shows that the secondary instability of the mixing layer 
is produced by spanwlse variations in the straining field of the primary 
vortex structures . 
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Chapter 1 


INIRODUCTION 


1 .1 Background 

Turbulent flows have been extensively studied, both experimentally 
and theoretically, for more than 50 years. The results of this effort 
still leave much that is unknown about the physics of turbulence, and 
our ability to predict its effects is still very limited. The theoreti- 
cal difficulties are primarily due to two things ; the nonlinearity of 
the governing equations and the large range of scales involved in high 
Reynolds number flows. The nonlinearity severely limits purely analyt- 
ical solutions. The wide range of scales cannot be handled numerically 
by present or foreseeable computers, and hence full 3-D, time-dependent, 
computational solutions can be obtained only at very low Reynolds num- 
ber . 

Large eddy simulation (LES) offers one approach to escape these 
difficulties. The LES approach is to mathematically distinguish "large" 
and "subgrid" scale components of a turbulent flow field. The equations 
for the "large" scales can be derived by smoothing or "filtering" the 
Navler-Stokes equations; but, due to the nonlinearity, "subgrid" scale 
terms appear, and these must be modeled. The merit of the LES approach 
(which is applicable to high Reynolds number flow) is due to the 
concentration of energy in the large scales or, correspondingly, to the 
fact that the modeled terms have a relatively small fraction of the 
energy and the large eddies do most of transport. 

LES is conceptually quite different from conventional phenomenolog- 
ical turbulence modeling, in which models for all scales of turbulence 
are required. See Reynolds (1976) for a review of those approaches. The 
limited success of modeling the entire turbulence field is understand- 
able when one considers the lack of universality in the large scale 
motions , since experimental evidence shows that large scale structures 
are very different in different flows. The work of Kline et al. (1959) 
on boundary layers and of Roshko (1976) on free shear flows illustrates 
this point. It is generally believed that the small scales tend to be 
universal in structure and hence should be much easier to model. 
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It Is also generally believed that large-scale structures contain 
most of the turbulent energy production and produce most of the impor- 
tant effects. The large scales are clearly dominant in transition of 
free shear layers (this can be proved analytically) and many researchers 
now believe this is also true far downstream (for example, see Browand 
and Troutt (1980) and Roshko (1976)). 

Large eddy simulation is thus an important tool for obtaining de- 
tailed Information about the most important turbulent motions. By 
examining the results of LES calculations, one can develop a deeper 
understanding of the physics of these flows. This insight should be 
useful in finding ways of controlling (to some degree) the effects of 
turbulence. Thus we believe that the most Important aspect of LES goes 
beyond prediction and into the domain of control; Liepmann (1979) has 
also pointed out the importance of the large structures in controlling 
turbulence. LES can also help in guiding the development of simpler 
phenomenological models by providing quantitative "data" for terms that 
must be modeled in these theories. 

The most Impressive LES work to date has been the turbulent channel 
flow simulation of Moin and Kim (1981) . This work clearly demonstrates 
the power of LES in providing both physical understanding and guidance 
for model development . 

A computational approach related to LES is Pull Turbulence Simula- 
tion (FTS), in which all of the scales of turbulent motion are computed 
in a 3-D, time-dependent calculation. With present computers, FTS is 
limited to low Reynolds numbers, when the large and small scales do not 
differ by more than one or two orders of magnitude. Nevertheless, FTS 
calculations are very useful, for they provide a testing ground for 
subgrid-scale models required in LES simulations, as well as physical 
insight and other "data" that are useful in developing phenomenological 
turbulence models. Examples of FTS calculation are Clark et al. (1977), 
Shlranl et al. (1981), Feiereisen et al. (1981), and Orszag and Patter- 
son (1972) and Rogallo (1981) 

The present work deals principally with LES computation of a time- 
developing mixing layer. This required development of a new numerical 
method. The results provide new insight into the physical processes 
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that occur during the transition of this layer from laminar to turbulent 
flow. 

1 .2 Related Work 

While the present work deals primarily with the LES approach, we 
shall now mention other computational approaches to gain perspective . 
In Reynolds' (1976) review, statistical methods are placed in a hier- 
archy ranging from algebraic to those in which five differential equa- 
tions are used to model the turbulence. Following a discussion of these 
statistical approaches, vortex methods will also be mentioned. 

Statistical approaches to treating turbulence generally employ a 
time average of the Navier-Stokes equations; occasionally other types of 
averaging are used. Time averaging the Navier-Stokes equations results 
in equations which govern the mean (time-averaged) velocity field. 
These equations contain higher-order terms — ‘time-averaged products of 
the fluctuating field. By averaging moments of the Navier-Stokes 
equations , governing equations for the time-averaged products of the 
fluctuating field are obtained. These equations contain averages of 
triple products of the fluctuating field. This process always results 
in higher-order terms — more unknowns than equations regardless of how 
many moments one takes . This is known as the closure problem of turbu- 
lence . 

To close the system of equations, a model or closure assumption 
which expresses the hlghest-order terms as some function of the lower- 
order terms is required. The lowest-order model, termed a zero equation 
model because no differential equations for the turbulence are employed, 
algebraicly relates the turbulent stresses (Reynolds stresses) to the 
mean flow quantities. Modeling at this level often produces satisfac- 
tory results for the mean flow, and calculations require only a few 
seconds of computer time . The success at this level is due to the 
considerable empirical content. 

The next level of modeling — a one-equation model — generally employs 
a differential equation for the turbulent kinetic energy and a 
prescribesd length scale or dissipation rate. This level is limited by 
the skill of the prescriber . 
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Two-equation models--some of which use equations for the turbulent 
kinetic energy and one which uses a length or dissipation scale — are 
very popular. They are still postdictive, as considerable empirical 
input is required to obtain good results; but they can be made to work 
well for given classes of flows . 

Full Reynolds stress models are used on occasion. At this level a 
differential equation is used for each component of the Reynolds stress 
tensor. While there is hope that such a model could be generally appli- 
cable, it is quite complicated. There are many empirical constants, and 
it is often hard to evaluate them. 

Statistical methods can "predict" the statistical behavior of the 
mixing layer. The present study, however, is aimed at better under- 
standing of the physics of the flow, and statistical approaches are not 
applicable. Approaches which are applicable are now presented. 

Ashurst (1979) has used a vortex tracking method to study the spa- 
tially developing mixing layer. In this approximation, the vorticity is 
discretized into two-dimensional line vortices which are periodically 
released at the origin of the mixing layer. The filiments are convected 
by use of the Blot-Savart law. A random perturbation of the initial 
release of the vortices causes the vortices to coalesce to form large- 
scale vortex structures, which pair. An Impressive color motion picture 
of this process was produced. The discretization approximation used is 
reasonable, though the 2-D approximation is quite limiting. The most 
ad-hoc assumption is the method of introducing artificial dissipation. 
This is necessary to prevent an unphysical growth in the kinetic 
energy. This approach is quite simple, fairly costly, and has an 
advantage over LES in terms of simplicity of the inflow and outflow 
boundary conditions and capability of treating a large region of the 
spatially developing mixing layer . 

This approach can be extended to three-dimensions, as shown by 
Leonard (1980) for a boundary layer. The full three-dimensional 
calculations are very costly, and we shall now discuss the vortex in 
cell method, which greatly reduces the cost. 

Couet and Leonard (1980) used the vortex in cell method to perform 
a three-dimensional simulation of the time-developing mixing layer . In 



this approach, three-dimensional vortex filaments are tracked. In place 
of the costly Blot-Savart method of tracking vortex filaments, the 
vortlcity is interpolated onto a fixed grid. The grid values of the 
vortlclty are used as a source for a Poisson equation for the vector 
potential whose curl gives the velocity field at the grid nodes. This 
velocity field is then Interpolated to the vortex filaments' positions, 
and the vortex filaments are convected to new locations. The cost is 
much lower than that for the Biot-Savart method. Couet and Leonard 
(1980) has developed a method of analytically continuing boundary 
conditions at infinity onto the edges of a finite grid, which eliminates 
image-flow problems (see Appendix D) . 

Couet's work produced good visualizations of the vortex formation, 
the existence of a secondary instability, and vortex pairing. He also 
studied the energy history of individual turbulence modes, but there is 
concern that the numerical methods used did not treat the conserved 
properties correctly. 

In an FTS, Metcalfe and Riley (1980) simulated a turbulent mixing 
layer. !l!heir full simulation required a 64"^ grid for adequate resolu- 
tion in a stream and spanwlse domain the same size as in the present 
work. They were restricted to a layer thickening by a factor of five or 
six, due to uniformly spaced grid in the gradient direction (the direc- 
tion of the mean velcoity gradient). Other than the treatment of the 
gradient direction, their numerics are very similar to those of the 
present work. Metcalfe and Riley's work demonstrates the capability of 
a full simulation to predict reasonably the mean velocity profile, 
growth rate, and turbulent energy of a time-developing mixing layer. 

Now we shall turn our attention to experimental work on mixing lay- 
ers. This is by no means a complete listing, but rather a few works on 
different aspects of the mixing layer. Reference to particular results 
relevant to the present work are given in Chapter 5. 

In a study of nonlinear eigenfunction interactions in transition, 
Miksad (1972, 1973) forced a laminar mixing layer at two forcing fre- 
quencies. Bradshaw (1966) studied the effect of initial conditions on 
the near field of a one-stream mixing layer by altering the geometry. 
Good sources of turbulence statistics in a mixing layer Include : Spencer 
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and Jones (1971), Wygnanski and Fiedler (1970), Patel (1973), and Cham- 
pagne et al. (1976). The as 3 nnptotic growth rate of the mixing layer is 
reviewed by Birch (1980). Visualization of the mixing layer Includes 
the work of Brown and Roshko (1974), Winant and Browand (1974), and 
Chandrsuda et al. (1978), 

1 .3 Motivation and Objectives 

The basic motivation for the present work is the desire to better 
understand the Importance of large scale structures in free shear flows 
through study of accurate, 3-D, time-dependent simulations of the for- 
mation of such structures. Most previous computational work on this 
problem was 2-D and often only the linearized equations were used. 
Recently, Metcalfe and Riley (1980) published one of the first three- 
dimensional simulations. To compute accurately the transition process 
accurately up to the point at which the energy reaches the experiment- 
ally measured levels for "fully developed" turbulence, the full non- 
linear equations must be used. Due to the fact, established analytic- 
ally, that transition or "near field” turbulence in the free shear layer 
is strongly dominated by a relatively small range of scales, this is a 
problem well suited to numerical simulation. 

To perform a simulation, we require knowledge of the flow field 
at the boundaries of the domain of interest. Often, this is a major 
stumbling block, as the boundary conditions at one or more boundaries of 
the computational region may be unknown. 

To Illustrate this point, consider the plane mixing layer. In the 
laboratory flow (see Fig. 1.3.1), two fluid streams moving at different 
velocities are initially separated by a splitter plate. At the end of 
the splitter plate the two streams form a free shear layer. This shear 
layer thickens with downstream distance. 

Experiments and theoretical analysis show that transition occurs 
through growth of a primary (Kelvin-Helmholtz) instability. This leads 
to the formation of spanwlse vortex structures. These in turn create a 
strong straining field between adjacent vortices, and in this region a 
secondary Instability leads to the formation of counter-rotating vorti- 
ces aligned with the strain. Pairing of the spanwise vortices is an 
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Important: mechanism in the formation of structures of larger scale. Far 
downstream! the mean velocity profile becomes self-similar, although 
quite different rates of growth have been found in different experi- 
ments . 

To simulate this flow, time-dependent inflow and outflow boundary 
conditions need to be imposed. Laboratory experiments never provide 
these in sufficient detail. Artificial conditions could be Imposed, but 
then there would be several grid points near both the inlet and the out- 
let where the computed solution would be unphysical. Hence a very long 
computational domain would be required; this approach is possible on the 
largest computers , but at very great cost . 

Tb avoid these problems with boundary conditions, we have chosen to 
treat a simpler problem, namely, the time-developing mixing layer (see 
Fig. 1.3.2). In the time-developing mixing layer, two semi-infinite 
streams of fluid travel in opposite directions with equal speed. The 
shear layer at the interface of these two streams thickens in time. 
Thiis flow is statistically homogeneous in planes parallel to the inter- 
face. Hence, periodic boundary conditions can be applied in the two 
coordinate directions in the plane of homogeneity. 

Because the layer is Immersed in an infinite region, some means of 
handling boundary conditions infinitely far above and below the layer is 
required. This could be done by applying free-stream conditions at a 
large but finite distance from the layer. However, this introduces 
"imaging errors” (Appendix D) . Instead, we shall use a method involving 
a non-uniform grid that extends to infinity above and below the layer. 
This requires the development of a new numerical method, which repre- 
sents one of the important contributions of this work. 

Let us consider the relationship between this idealized flow and 
the laboratory flow. The laboratory flow is characterized by a param- 
eter A, defined by 


X - (Uj - n2>/(Uj + Uj) 

where Uj^ and U 2 are the speeds of the two streams . As A goes 
to zero, all streamwlse mean gradients become smaller, and the space- 
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developing layer becomes more like the time-developing one. Unfortu- 
nately, when X is small, the early part of the flow is dominated by 
the wake of the splitter plate and is not really a mixing layer flow. 
Hence, we cannot make a direct comparison of the simulation with experi- 
ments. Nevertheless, we believe that the basic physical processes are 
the same in both flows , and hence can be examined by study of the time- 
developing flow. The most significant exception is that, in the labora- 
tory flow, the downstream flow can affect the upstream flow by creating 
pressure fluctuations. This allows "phase locking," a situation in 
which the flow remains in a particular configuration for a long time. 

Simulation of the time-developing flow requires initial conditions. 
The effect of different initial conditions in the time-developing flow 
is probably similar to the effect of different upstream conditions in 
the laboratory flow. The sensitivity of the near field of the mixing 
layer to initial conditions was examined in the laboratory by Bradshaw 
(1966) . He found the near field to be quite sensitive to changes in the 
Initial disturbance field. In this laboratory study, the initial dis- 
turbance field was altered by geometry changes affecting the boundary 
layer of the one-stream mixing layer. 

In the computer simulation we can specify any initial field that is 
consistent with a few simple constraints , and hence we can examine the 
effects of changing specific features of the initial disturbance. All 
computational cases will begin with the same mean field; for each case, 
a specific disturbance field will be added to complete the initial 
condition. Changes in the initial disturbance field to be considered 
Include: amplitude, spectrimi shape, and the set of random phases of the 
modes . 

After a sufficient period of development, we expect that the mean 
(phase- averaged) velocity profile will become self-similar. The rate at 
which this occurs , and the degree of self-similarity in other turbulence 
parameters of the time-developing flow should provide insight into the 
self-similarity (or lack thereof) in the laboratory flow. 
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Summarizing, the objectives of this work are thus: 

• To develop an accurate method of differentiation and inte- 
gration on a non-uniform grid, suitable for simulating a 
time-developing mixing layer. 

• To employ the above method in a study of simulations of the 
transition and early turbulence domains of a time-developing 
mixing layer. 

1 .4 Summary of Present Work 

A new and very accurate numerical differencing and Integrating 
scheme for infinite domains is presented. It is based on the use of 
Fourier expansions and takes advantage of the computational efficiency 
of the fast Fourier transform. The new method is applicable to more 
general boundary conditions than the standard Fourier method, due to the 
use of mapping functions. (The simplest boundary conditions to imple- 
ment are periodicity , or zero , or zero-derivative conditions , or combi- 
nations thereof.) However, the allowed mapping functions are restricted 
for reasons of efficiency and accuracy. For more detail, see Chapter 3. 

Two particular mapping schemes, both for doubly infinite domains, 
are implemented. One is chosen to handle jet-type flows, while the 
other is designed for the mixing layer. Both schemes are applied to 
linear test equations having known analytical solutions. The new scheme 
is shown to have errors as much as six orders of magnitude smaller than 
common finite-difference schemes for equal numbers of mesh points. 

The new scheme designed for jet-type flows is used to demonstrate 
the Influence of finite-domain boundary conditions. The evaluation was 
made by use of an array of two-dimensional vortex structures , all with 
the same sign of vortlclty. If this array of vortices is computed in a 
finite domain with no-stress boundary conditions applied to a surface 
parallel to the array, but at a finite distance, image flows are im- 
plied. The nearest image flow above is the mirror image (Imaged about 
the no-stress surface) of the initial vortex array. This image vortex 
array and the true vortex array form a near-field jet-type flow, and 
their beliavlor in time was computed using the new infinite-domain 
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scheme. Comparing a base computation with one In which the Image array 
was shifted In the streamwlse direction (relative to the true array) 
showed potentially strong coupling to Image flows , depending on vertical 
spacing. For details, see Appendix D. 

Using the new Infinite- domain scheme, a 3-D, time-dependent, large 
eddy simulation study of transition and early turbulence In a time- 
developing mixing layer was undertaken. The primary focus of this study 
concerns the effect of the Initial disturbance field on turbulence 
development. Effects due to filtering and modeling are also examined. 

To sort out the effects of the Initial disturbance field, the same 
laminar, mean-velocity profile Is used as the Initial mean field In all 
cases. Tb this mean velocity field, an Initial divergence-free distur- 
bance field Is added. We use nine cases Involving seven different Ini- 
tial disturbance fields. These seven cases examine the Influence of the 
disturbance amplitude, spectrum shape, and random phase sets on the 
resulting early turbulence. 

The computations provide the mean velocity profile, the momentum 
thickness, the turbulent kinetic energy, the Reynolds stress tensor, the 
Reynolds stress anisotropy tensor, and particle tracking pictures. 
Examination of these results provides better understanding of the mixing 
layer. The central observations from these computations are as follows : 

• All cases display Immediate self-similarity in the mean vel- 
ocity profile. 

• The momentum- thickness growth rate is strongly Influenced by 
the initial disturbance-spectrum shape. 

• Interesting oscillatory behavior is observed in the kinetic 
energy profile width for the small-amplitude initial distur- 
bance cases. TMs oscillatory behavior is not present in 
high- amplitude initial disturbance cases. 

• The anisotropy tensor proved to be the most sensitive mea- 
sure of self-similarity. Even changing the random phase 
distribution in the initial disturbance field produced 
enormous differences in the evolution of the anisotropy 
tensor . 

Probably the most significant aspect of the study was revealed in 
the particle-track pictures. While large coherent structures readily 
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appeared (some similar to those of Winant and Browand (1974) and others 
like Chandrsuda et al. (1978) j the mechanism for producing the secondary 
Instability has been identified. It seems to be the result of spanwlse 
variations in the strength or position of the primary vortex structures , 
which give rise to spanwlse variations in the straining field stagnation 
line. This causes spanwlse vortlclty to be tilted towards the stream- 
wise direction as the vorticity is rolled up by the primary vortices , 
and this process leads to the formation of pairs of counter-rotating 
vortices aligned with the straining field . 

Changing the set of random phases of the initial disturbance field 
produced significantly different results, both in the statistical and 
structural characteristics of the mixing layer. This is a consequence 
of the small sample of large eddies in any given calculation. The im- 
plications of this are still very significant, however, as any given 
experimental apparatus is likely to produce a given type of large eddy 
structure, which locks on. This is due to the likelihood of a fixed 
type of initial disturbance being present and also due to pressure feed- 
back effects. However, different experimental apparatus is likely to 
lock onto different large eddy types. There are two and perhaps more 
different large eddy patterns which are possible . Thus , large eddy 
variation may be small in a given experiment, but significant variations 
may occur from experiment to experiment . 

TWo cases were run to examine the effects of filtering and subgrid 
turbulence modeling. We found that filtering delays the onset of non- 
linear effects and gives us less than the total picture. However, it 
considerably extends the length of time over which the computation is 
meaningful. The subgrid-scale model was shown to have very little 
influence on the calculations . 
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Chapter 2 


PROBLEM FORMULATION 


2 .1 Governing Equations 

We shall restrict ourselves to considering the motion of an Incom- 
pressible Newtonian fluid. The motion of such a fluid is governed by 
the Mavler-Stokes equations. A common form of these equations Is: 


3Uj 


1 3p 
P 


VV^U, 


(2.1.1a) 


For computational conservation of important properties the equations 
are written In the form : 


«Ui jp 

“j -5^" "j ^ 

The modified pressure is P = (p/p) + (Uj^Uj^/2) . 
servatlon of mass of an Incompressible fluid Is : 


vV^u^ (2.1.1b) 

The equation for con- 


9u. 

* 0 (2.1 .2a) 

For computational convenience we take divergence of (2.1.1b) and 
enforce (2.1.2a) to get a Poisson equation for the modified pressure: 


9^ 

v^p 


3u. 3u. 3u. 3u . 3^u . 

J 

1 


(2.1.1b) 


Thus, given the velocity of a fluid field at some time t^, we may 
solve (2.1.2b) for the modified pressure and then find the time rate Of 
change of the velocity field from (2.1.1a). From this, we can find the 
velocity at the next time step. 


This form allows many differencing schemes to conserve mass , mo- 
mentum, and energy (as shown by Mansour et al., 1977). 
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2.2 Boundary Conditions 


The boundary condition Imposed on the time-developing flow In di- 
rections of statistical homogeneity Is that all variables be periodic : 


Uj^Cx.y.z) 
u^Cx.y ,z) 


Uj.(jd-Lj^,y ,z) 
Ui(x,y+L2,z) 


(2.2.1a) 

(2.2.1b) 


In the direction of Inhomogeneity, we can require that the flow become 
the unperturbed free stream far from the shear layer, or, alternatively, 
we can impose the no-stress conditions : 


■^(x,y,z)|^^ = 0 = (x.y.z)| 2 =±eo » w(x,y,±~) = 0 (2.2.1c) 


The no-stress condition Is advantageous numerically and Is employed In 
the present work. 


2 .3 Initial Conditions 

The Initial velocity field used consists of a laminar field plus 
small perturbations. The time-developing laminar mixing layer has 
gradients only in the z direction, and only the u component of the 
velocity field is nonzero. This layer is thus governed by the diffusion 
equation : 


3u 

It “ 



(2.3.1) 


With the no-stress boundary condition and Initial condition u(z,0) = 
u^ for z > 0 and u(z,0) = -u^ for z < 0, Eq . (2.3.1) has the 
solution 


u(x,y,z,t) = u erf(z/Avt) 
o 


(2.3.2) 


The slope thickness, 6^, of a free shear layer is defined as the 
velocity difference across the layer divided by the maximum gradient of 
the layer . For the laminar layer ; 
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(2.3.3) 


6 = v^4irvt 

u 

For our study of transition, we shall start our calculation at a time at 
which Reg = =60. 

The disturbance field is constructed in the following way. First, 
a periodic, divergence-free, homogeneous, isotropic field on a 16 x 16 
X 16 grid is constructed using the routine written by Kwak et al. 
(1975) . The velocities on five planes of this field are assigned to the 
central five planes of the grid used in the present calculation (which 
is non-uniform and anisotropic) . This field is then smoothed in the 
gradient direction by a Gaussian filter; the result is smooth, but not 
divergence-free. We obtain a divergence-free field by taking the curl 
of this field. The divergence-free velocity field is added to the 
error-function profile to give the complete initial field. Further 
details of this process are given in Chapters 3 and 5 . 

2 .4 The Computational Domain 

The choice of a streamwlse and spanwlse box length (Lj^ and I^, 
respectively) is critical. Michalke (1965) and Betchov and Crimlnale 
(1967) studied the stability characteristics of a mixing layer with a 
hyperbolic tangent profile. While our error function profile is not 
identical to theirs, it is sufficiently close that we can use their 
results as a guide. Betchov and Crimlnale considered the linearized 
Navler-Stokes equations (valid for small amplitude disturbances) and 
searched for elgensolutlons growing in time. They found, even in the 
limit of infinite Reynolds number, that there is a minimum wavelength 
in the streamwlse direction for which the disturbance is amplified 

by the mean shear; any disturbance with a shorter wavelength will decay. 
In non-dimensional terms at infinite Reynolds number, the shortest 
amplified wave has a wavenumber 

a = = 1.0 (2.4.1) 

c A, 

Ic 

According to Betchov and Crimlnale, the most amplified disturbance wave 
length at infinite Reynolds number corresponds to wavenumber = 0.43. 
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The effect of finite Reynolds number is to decrease the values of 
and at the Reynolds number we shall use, = 0.34. Using these 

considerations as a guide, we set the computational domain such that the 
smallest a supported in the calculation was 0.17. This corresponds to 
the longest wave allowed — the one whose wavelength is the computational 
box length. Since we are using Fourier methods, we have adequate reso- 
lution and a large enough domain to study the most amplified waves. At 
our initial Reynolds number (Reg = 60) and with a grid spacing in the 
span direction equal to the strearawlse grid, there are 14 amplified 
modes in our discrete approximation (this includes 3-D modes) . 

2 .5 Filtering 

In treating a turbulent flow numerically, we may have more scales 
of motion than any computer can handle; this depends somewhat on the 
flow and Bteynolds number. If this is the case, we are forced to filter 
the flow field in a way which leaves a range of scales that can be com- 
puted. We follow Leonard (1974), who first formalized this approach. 
Since we know from experiments that the largest scales of motion are the 
most energetic ones and are responsible for most of the transport, we 
shall truncate the small scales . We shall symbolize the large-scale 
field by u. For a homogeneous flow, we define it as: 

u(x) = ^ G(3^x’ ) u(x') dx' (2.5.1) 

(integration over all space) . The small-scale or subgrid (SGS) field is 
simply the difference between the full field and the large-scale field. 

u* = u - u (2 .3 .2) 


Note that, since IT u, TT*” 0. 

The choice of a difference kernel for the filter function has two 
assets. Most Importantly, such a filter commutes with differentiation 
operators . This means that both the large-scale field u and the 
subgrid scale (SGS) field u' will separately satisfy the continuity 
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equation. The second asset of a difference kernel for the filter re- 
sults from the convolution theorem. Discrete Fourier methods can be 
used In our problem, and the convolution theorem allows a very fast way 
of computing u by means of the Fourier transform. The transform of 
u Is : 

A 

A A 

u(k) ■= G(k) u(k) (2.5.3) 


and u Is obtained by Inverting (2.5.3). 

Our choice for G(x-x') Is a Gaussian (see Mansour et al., 1977): 


G(x-x') = exp 


(2.5.4) 


with = 2h^, where h^^ Is the grid spacing In the 1*-*^ direction, 
and n Is the number of directions In which we elect to filter. It 
should be noted that the discrete transform of (2.5.4) Is used In 
(2.5.3). This Is slightly different from the continuous transform. 


We are using a non-uniform grid scheme In the z or x^ direc- 
tion. Moln et al. (1978) showed that It Is difficult to define a filter 
for this direction. One can Introduce an approximation, but the large 
and SGS flow fields would no longer be divergence-free. Thus we do not 
filter at all In the x^ direction. If we apply the filter to (2.1.2a) 
In only the Xj^ and X2 directions (n = 2) , we get 



(2.5.5) 


Recalling u^^ = u^ + uj^, Eqs. (2.1.2a) and (2.5.5) Imply that the 
subgrld-scale field Is divergence-free as well. This demonstrates the 
desirability of a linear filtering operator which commutes with differ- 
entiation. 


Application of the filter to (2.1.1b) gives 


3 - 
It ^1 


- ®“i 

“jHj 



3P „2— 
+ vV Uj^ 



(2.5.6) 


where 
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- \k 


jLi 

3 


TT P 1 / \ 

p = L+-^+^(u^U^) 


\j = u’Uj + u^u' + TJfSy 


Note that contains subgrid-scale terms and therefore must be mod- 

eled. 

Taking the divergence of Eq. (2.5.6) and enforcing Eq . (2.5.5) 
gives a Poisson equation for the filtered modified pressure ; 


2- ^“i 

V^P = — 1 — i 




a"- - 3 3 


ij 


(2.5.7) 


Since filtering is a linear operator and all boundary conditions 
are linear , we arrive at the following boundary conditions on the large- 
scale field: 


uj^(x,y,z) = uj|^(x + Lj^.y.z) 


(2.5.8) 


Uj^(x,y,z) = Uj^(x,y + L2,z) 


(2.5.9) 


|| (x.y.z) 


3v 


z=±“ 


0 = ^ (x.y.z) 


; w(x,y.±«») = 0 (2.5.10) 


z=±« 


2.6 Subgrid Modeling 

Since contains subgrid-scale terms, it must be modeled. The 

history of how we came to the model eventually used is of some Interest . 

In a preliminary phase of this work, we explored the fully tur- 
bulent mixing layer. We solved the vorticity equations using the vor- 
ticity model developed by Mansour et al. (1978). We also used the 
primitive equations with the following model, due to Smagorinsky (1963) : 


iJ 


- . 

T ij 
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( 2 . 6 . 1 ) 


hi - 

In (2.6.1), A Is the filter width and C„ Is a constant, approxl- 
mately 0.2. Runs with the same Initial conditions — one using the prim- 
itive equations and Smagorlnsky's model and another using the vortlclty 
equations and Mansour* s model — had turbulent statistics which were 
essentially Identical. In both runs, however, the mixing layer thick- 
ened much faster than the experimental layer . In an attempt to under- 
stand why, we then used Smagorlnsky's model In a calculation of an 
Initially laminar mixing layer with no disturbance added , and found that 
the layer grew between two and three times as fast as the experimental 
turbulent layer . This erroneous behavior provided the clues needed to 
modify the model, and led to an Improved model that we used In the 
transition studies. 


One reason for the Improper behavior In laminar flow arises because 
the model "turns on" too quickly. The models are supposed to account 
for the effects of unresolved turbulence , but the eddy viscosity con- 
tains a significant contribution from the mean field. Hence, In a 
laminar flow with no turbulence, these models Incorrectly provide eddy 
viscosity and hence subgrid stresses . A model that allows the subgrid 
stresses to build up slowly with the turbulence field can be made by 
redefining v,j, as 




f 

4 


2(S 


ij 


< S,, »(S,, - < S,, » 


ij 


ij 


ij 


-] 


1/2 


( 2 . 6 . 2 ) 


where the < > denotes an average over a plane of constant z . 

All model calculations reported here use given by (2.6.2). 
Another modeling concept Is presented In Appendix C; work by Bardina et 
al. (1980) suggests that this new model Is very promising. 

The problem with too-rapld growth of the turbulent layer was also 
related to the grid layout. The stability considerations discussed In 
Section 2.4 require a grid spacing that Is very large In the streamwlse 
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and spanwise directions (the only directions in which we filter); in 
fact, the spacing is large compared to the thickness of the layer, and 
hence the filter width was much larger than the layer thickness. Phe- 
nomenological models typically use length scales for the large eddies 
that are about one-tenth the shear layer thickness; this is far shorter 
than our filter width. Hence, we limited our length scale A to a max- 
Imtun of one-tenth of the shear layer thickness A^ = min(2h^ ,6^/10) . 

2.7 Summary 

In summary, in the mixing layer transition study we solve the equa- 
tion for the filtered field u (2.5.5 and 2.5.6) using the subgrid 
model (2.6.1) with from (2.6.2). The filter (2.5.4) was used, and 
its width Aj. was 2hj^. The boundary conditions (2.5.8-10) were ap- 
plied. The initial velocity field consisted of a laminar mean field 
(2.3.2) with Reg = 60, plus a divergence-free random perturbation (see 
Chs . 3 and 5). The computational domain was chosen to allow a number of 
the amplified modes of the laminar Instability to appear in the solu- 
tion, including the most rapidly growing mode. 

At this point, the global problem formulation is complete. We now 
proceed to the details of the numerical method used. 
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Chapter 3 


NUMERICAL METHODS 

3 ..1 Preview 

Some familiar and some new numerical methods are used In this work. 
To orient the reader, we shall briefly preview the methods to be dis- 
cussed in this chapter. 

Our choices of numerical methods were guided by the objectives of 
this study. Desirable methods preserve as much of the physics as pos- 
sible, and this requires accurate numerical representation of the spa- 
tial derivatives. Fourier methods provide the most accurate differ- 
entiation for a given number of grid points; hence we used them wherever 
possible . 

Since periodic boundary conditions are applied in two directions, 

and X 2 ) , we can use the standard Fourier scheme (described in the 
next section) to treat spatial derivatives in those directions. The 
gradient (x^) direction requires special consideration; we were con- 
cerned about the influence of image flows , which are discussed in Appen- 
dix D. lb avoid imaging problems, we chose to compute the solution over 
infinite x^ . For reasons of accuracy we wanted a discrete orthogonal 
function expansion method which could treat an infinite region. Since 
no existing method was known, we developed a mapping scheme which re- 
tains the efficiency of the fast Fourier transfom (FFT) and is ideally 
suited to our problem as well as several others. This method is de- 
scribed in Sections 3.3 and 3.4. 

New measures of spatial resolution and statistical validity in 
Fourier methods were developed (Section 3.5) and used to quantify the 
accuracy and statistical validity of the calculations . These are de- 
scribed in Section 3.5. 

To aid the understanding of the physics , computational flow visual- 
ization is used. The visualization is achieved by tracking the inter- 
sections of a freely deforming grid. While this is similar to particle 
tracking, it differs in that the identity of the individual particles is 
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retained. Interpolation of the velocity field is required for tracking 
the intersections, and this is presented in Section 3.6. 

Recall that the initial field consists of a laminar mean field 
plus a random disturbance. The method of constructing the initial 
divergence-free disturbance velocity field is described in Section 3.7. 

Three time-advance schemes were used in parts of this study. The 
transition problem has a disturbance velocity field which grows by 
orders of magnitude during the simulation. For this reason, a time- 
advance scheme which allows the time step size to change continually 
will be advantageous. To maintain accuracy, we want an explicit scheme 
of fairly high order. The fourth-order Runge-Kutta method (Appendix 
B.2) is ideal. It satisfies the above conditions and has excellent 
stability characteristics as well. All of our numerical studies of 
transition used this method for time-advance of the velocity field. 

To advance the particle positions in time , we require less accu- 
racy, since this is used for visualization purposes only. We still 
require a scheme which allows continual adjustment of the time step and 
has moderate accuracy. The second-order Runge-Kutta scheme (Appendix 
B.3) was chosen for this problem. 

In addition, some computational experiments discussed in Appendix D 
were carried out using the second-order Adams-Bashforth method for time 
advance (Appendix B.l). 

3 .2 Standard Fourier Methods 

The discrete Fourier transform is defined for any number of grid 
points; however, efficient algorithms usually limit the number of points 
to particular Integers. In the most popular routines due to Cooley and 
Tukey (1965), this number must be a power of two; the Winograd (1976) 
method allows other numbers to be used , but has not yet seen extensive 
application. Thus, we consider only variations of the Cooley-Tukey 
algorithm. 

If we have a function defined on a set of uniformly spaced grid 
points, say Xj^ = (n-l)/N, the discrete Fourier transform can be de- 
fined by : 
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A ik2trx 

\ - E 

n“l 

We shall assume N = 2™. There is a great deal of freedom in the choice 
of wavemumbers k used in Eq. (3.2.1). (Tills formulation gives inte- 

ger wavenumbers) . So long as we are interested only in representing the 
function on the grid points, the choice (within the allowed bounds) is 
Irrelevant. However, when we use the fast Fourier transform (FFT) as a 
means of obtaining derivatives we are regarding Eq. (3.2.1) as an inter- 
polation, and it is Important to choose the wavenumbers which give the 
smoothest interpolation possible. For an even number of grid points, 
there are two equally good choices : 


or 
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Either of these choices means that the wavenumber |k| = N/2 is 
represented only by a single waveform, whereas all other wavenumbers 
have two waveforms, l.e., ± k. As a result, we have Incomplete informa- 
tion about the highest wavenumber component, in fact, we know neither 
its phase nor its amplitude. Consequently, we cannot differentiate it. 
Most workers set its derivative equal to zero to avoid the problem. 
(This problem does not occur with odd-point transforms.) The derivative 
is thus : 


where 
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(3.2.2) 
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3 .3 Fourier Transforms on Non-Uniform Grids 

The usual method of dealing with an Infinite region Involves map- 
ping It onto a finite region. Mappings are also commonly used to modify 
the geometry such that the function Is smoother In the transformed coor- 
dinate system. Mappings Invariably complicate the equations to be 
solved , but they offer the Important advantage that numerical methods 
are both more easily applied and more accurate In the transformed coor- 
dinate system. These advantages almost always outweigh the disadvan- 
tages, and coordinate transformations are a standard part of numerical 
methods today. Indeed, the development of better mappings Is a major 
field of research. 

In describing the new method, we shall restrict our attention to 
one-dlmenslonal problems . Suppose that z Is the physical coordinate 
and we Introduce the computational coordinate t, by means of the 
mapping 

z = h(C) (3.3.1) 

The derivatives In the two coordinate systems are related via the chain 
rule ; 


^ 1 If 

dz ~ d? dz TT” d? 


(3.3.2) 


Now, suppose that we choose to represent the function In the trans- 
formed coordinate system In terms of Its values on a uniformly spaced 
grid Cj = jA?. This function can be represented In terms of Its 

Fourier transform In the manner described In the previous section. 
Thus , we can write : 


N/2-1 

K 

k=-N/2 


lk2nC . 


(3.3.3) 


We could use this Fourier transform to compute the derivative 
df/dC that appears In Eq . (3.2.2). One could then substitute the 

result Into Eq. (3.3.3) and compute df/dz. The difficulties with this 
approach are (1) that. In general, the result contains a considerable 
aliasing or truncation error, and (2) that the resulting operator cannot 
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be Inverted when the number of points used Is even. Thus we must seek 
further improvements . 

To look into this more deeply, we begin by noting that l/h*(?), 
which appears in Eq. (3.3.2), can itself be represented as a discrete 
Fourier series similar to the one in Eq . (3.3.3). In general, N terms 
will be needed and, when this series is multiplied by the one represent- 
ing df/dC, the result will contain 2N wavenumbers, -N, ... , N-1 . 
Truncatiiog the result to N terms produces the truncation error alluded 
to above. 


The problem can be avoided by restricting the allowed mapping func- 
tions to those which contain only a few Fourier modes with small wave- 
numbers . Thus , we can write : 

m 


1 _ " ik2iic 

h '(0 Z-/ \ ® 

k=-m 


m « 


N 


(3.3.4) 


When this is substituted into Eq . (3.3.2), the result contains N + 2m 
wavenumbers . By making m small and truncating the modes whose wave- 
numbers are less than -N/2 +1 or greater than N/2 - 1, we produce a 
small, acceptable amount of truncation error; doing this accurately re- 
quires that the multiplication of 1/h' and df/d? be carried out in 
Fourier space. It is possible to take the product in configuration 
space; however, the result will be aliased and will populate the ± n/ 2 
modes an<l thereby make the Inverse of the differentiation operator sin- 
gular. Defining the derivative via the truncated transform- space prod- 
uct allows us to construct an Integral operator which is the "exact 
inverse" of the differentiation operator. By "exact inverse" we mean 
that the derivative of the Integral of f is exactly f. Note that 
these are alias-free operators . 


We now apply these ideas to two mappings suitable for free shear 
layer problems in fluid mechanics. The problems of Interest are best 
treated in infinite domains . For example , in the computation of plane 
jet flow the region of Interest is doubly infinite in the gradient 


See ^pendix D for a demonstration of the danger of using a finite 
domain to study vortex pairing . 
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direction, the boundary conditions are Identical at * and we would 
like the grid points to cluster near the origin. The cotangent is a 
suitable mapping function for this purpose , 1 .e . , 

z = h(5) = -a cotan(TT?) ; -“<z<“ (3.3.5) 


This gives the metric 


1 

h' 


^ sin^Cir?) 


1 

lia 


1 - 




(3.3.6) 


Recalling Eq . (3.3.4), we see that m = 1 and the Fourier coeffi- 
cients of the metric are : 




“-1 


1 

4na ’ 


1 

Ua 


(3.3.7) 


We thus have grid clustering near the origin and a minimiua of truncation 
error, since m = 1. An estimate of the error in the derivative of 
f(z) is given by 

This error will be small if the Fourier series for f(z) converges 

rapidly. The mapping above is applied to a vortex-pairing problem In 
Appendix D. It is also useful in treating time-developing plane jet 
flow. 

The problem of Interest in this work is the time-developing mixing 
layer. The main difference between this case and the previous one is 
that the boundary conditions are no-stress rather than periodic . A var- 
iation of the mapping given in Eq. (3.3.5) is suitable for this case, 
namely , 

z = h(5) = - a cotan(2iT5) , -“^z^“ (3.3.9) 

As indicated, the domain 0 1/2 is the image of the physical 

region - ^ z » under this mapping . However , the boundary condi- 

tions are such that the problem is not periodic in this domain. We 
shall therefore let ? range from zero to unity in the computation. 

The data are made periodic in this domain by defining the solution for 
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1/2 ^ 1 by reflection about C = 1/2. This means that In z-space 
we are considering two Rlemann sheets, l.e., we are essentially doing 
the problem twice. The difference between this method and the previous 
one is equivalent to using a Fourier sine or cosine expansion in place 
of an exponential Fourier transform in a finite domain; in either case, 
twice the work Is required. The choice of sine or cosine transform 
depends on the nature of the function being expanded. The metric re- 
sulting from this mapping Is : 


1 

h*- 



sln^(2ir?) 


1 


• ^^12(2x0 ^ ^-12(2vOj-| 

1 2 


(3.3.10) 


From Eq. (3.3.4), we now have m = 2, and the Fourier coefficients 
of the metric are : 




a 

o 


(3-3.11) 


The truncation error Is again small; an estimate for the maximum 
error in the derivative of f in this case is : 


e 


2Tf|a2 I 




(3.3.12) 


Equations (3.3.5) and (3.3.9) are mappings of an infinite physical 
domain onto a finite domain with grid points clustered near the 
origin. Equation (3.3.5) Is suitable for periodic boundary conditions, 
Eq . (3.3.9) Is suitable for no-stress conditions, and was used in our 
free shear layer simulations. 


3 .4 Derivative and Integral Operators 


It is Important In numerical analysis to use Integral and deriva- 
tive operators that are exact Inverses of one another , 1 .e . , that are 
"consistent." With the results of the last section, we can define a 

consistent set of derivative and integral operators. Using the mapping 
(3.3.9) and recalling Eqs . (3.3.2) and (3.3.4), we have 


dz 


1 

4ira 


2 - (e 


12(2ttc.) -l2(2TfC.) 

^ + e 
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N/2-1 
k=-N/2+l 


ik2Tr? 
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Truncating the above expression to include only the wavenumbers k 
between -N/2 + 1 and N/2 - 1 gives the definition of the first de- 
rivative operator : 


6z 


1_ 

2a 


N/2-1 


r 

k=-N/2+l 


i(k-2) ; 

2 ^k-2 


i(k+2) 

2 



i2uk?j 


(3.4.1) 


The prime on the summation indicates that any term whose subscript has 
magnitude greater than (N/2 - 1) is zero. 


The second derivative operator is obtained by a second application 
of the first derivative operator , giving ; 

N/2-1 




6z 


N/2-i ( 

, - k E' kK- 

'j k=-N/2+l ( 


i(k-2) 


i(k+2) 


"k-2 


^k+2_ 


i(k-2) 

4a 


i(k+2) 
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k-4 


ik 

T~ 


\ 

'k_ 


i(k-Hi) ; 

~2 k+4 


(3.4.2) 


i2irkc j 


We also a need a fast and accurate Poisson solver for infinite 
domains . The Poisson equation 

V2p = Q (3.4.3) 


may be solved by use of the three-dimensional Fourier transform. For 
the standard case of a function periodic in all three directions and a 
grid uniformly spaced in each direction, we may find the Fourier coeffi- 

A A 

dent P(k) by dividing the corresponding Fourier coefficient Q(k) by 
the negative square of the magnitude of the wave vector J^. One then 
inverts the Fourier transform to get P Itself. This is a very effici- 
ent and accurate method of solution. If we use a non-uniform grid in 
one direction (but uniformly spaced in the other two directions) , the 
solution procedure is only slightly more complicated. For illustrative 
purposes, we use the mapping (3.3.9) which leads to the second deriva- 
tive operator given by (3.4.2). Equation (3.4.3) may be solved by con- 
sidering the linear algebraic system : 


27 



The non-zero elements of A are : 


= _ (k-2)(k-4) 

.. k(k-2) (k-2)^ 

" 8a^ 8a^ 

(3.4.5) 

4a 16a 16a 

- k(k+2) (k+2)^ 

■■ g ^2 + 3^2 

„ « (k+2)(k-f4 ) 

k ,k+4 7~1 

16a 

In this case , kj^ and k£ are the wavenumbers in the uniform 
grid directions, while k is the C-wavenumber in the non-uniform direc- 
tion. Also, in (3.4.5) any factor whose magnitude is greater than (H/2 
- 1) is set to zero. (The problem for uniform grids in all three 
directions can be viewed as solving (3.4.4) when A is a diagonal 
matrix.) The solution to pentadlagonal system (3.4.4) can be quickly 
and easily solved on the computer. Note that, because of the differen- 
tiation/integration consistency, this formulation results in a Laplaclan 
identical to the divergence of the gradient operator. This condition is 
necessary to maintain a divergence-free velocity field, as Kwak et al. 
(1975) have pointed out. 



3 .5 Resolution and Statistical Validity 

This section will present measures that will help assess the valid- 
ity of the calculations. We shall use discrete orthogonal function esc- 
panslons in each spatial direction. This is the most accurate approach 
available. However, the nonlinear terms in the governing equations make 
the issue of accuracy difficult to deal with. The product of two vari- 
ables that can be represented using N modes requires 2N modes for 
complete accuracy. Since we take the product in configuration space, 
the result is aliased. Aliasing; is the pollution of the low wavenumbers 



due to high wavenumber Information masquerading as low wavenumber infor- 
mation when the grid is too coarse. 

To obtain a measure of the accuracy of a configuration space 
product, we first define the energy spectrxim: 


E(k^ ,k 2 ,z) 


1 “l"2 


The asterisk denotes the complex conjugate. In (3.5.1) we have Fourier 
transformed only in the plane of homogeneity, since the grid is non- 
uniform in z. Consider the central plane (z = 0) and note that the 
highest fully resolved mode in the i*-^ direction is the Nj^/2 - 1 
mode. In the product of u^ with itself, if all non-zero modes have a 
wavenumber magnitude |kj^| satisfying 


Ik^l < Nj^/4 (3.5.2) 

then the product is also fully resolved with no aliasing on a grid with 
points. The fraction of the flow computed with full resolution is 
therefore 


n. = 


- Z’ZT /L 

k. k- / k- k- 


0 ) 


The double prime on a summation indicates that Ikj^l < N^/4 . 
further the alias-free fraction of the energy is given by 


(3.5.3) 


We note 
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Ml 


E(k^ ,k2 ,0) 


k. 


k. 


/E£«- 


2^ ,k2 ,0) 


(3.5.4) 


The triple prime on the summation indicates the sum over all |kj^| < 
N^/4. One further comment is appropriate, namely, the above are sums 
over squares in k space. 


We now take up the issue of sampling. In a computation of tran- 
sition in a mixing layer, large, coherent structures are formed. Not 
all are identical , and thus we desire a measure of how many large 
structures we capture. We define an energy-weighted measure of the 
sample size : 
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(3.5.5) 


X/ max( Ikj^l ,1) max(|k 2 l,l) E(kj^,k 2 , 


0 ) 


Tci k/i 


E(k^ ,k^ ,0) 


The sample size tells us liow many effective full Fourier waves the 
total energy has in the computational domain. It is essentially the 
ratio of the area of the computational domain to the product of length 
scales in two coordinate directions; the definition of the length scales 
is implicit in Eq . (3.5.5). 


For a given number of grid points , one may choose a fine mesh so 
tlliat and are very near one, but the sample Tig, will be 
small. Conversely, with a coarse grid rig is large, but we sacrifice 
accuracy. In other words, we may compute the behavior of one large eddy 
accurately, or compute several large eddies crudely. The choice depends 
on one's objectives. Values for n^, ri^p, and Tig are presented with 
the computational results in Chapter 5. 


3 .6 Interpolation 

To gain further insight into the physical nature of transition in 
the mixing layer, we shall track particles using computer graphics. 
Aifter coiQputing the flow field, we know the velocity field at the grid 
points, and we must Interpolate to find the velocity at the location of 
each particle. For maximum accuracy and consistency, we should use 
Fourier interpolation. However, the cost of Fourier interpolation is 
excessive; for this purpose. We chose to use a three-point interpolation 
scheme (three points in each direction, or a 27-point box). In the uni- 
form grid directions , the function is represented as 


f(x) = 


C + C, X 

o 1 


+ CjX 




(3.6.1) 


The origin Is always taken to be the nearest grid point (which we call 
Xj.^) and the distance between grid points is taken to be unity. At the 
nearest grid point the value of the function is and we have 
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(3.6.2) 
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The interpolation In the z direction (non-uniform grid) is simi- 
lar, but more complicated. It is based on a Taylor expansion, and the 
first and second derivatives are evaluated using the formulas given in 
Appendices A.l and A. 3. Also, rather than using the nearest grid point 
as the origin in the z direction j we use the grid point for which 
h'(z^)|z - Zji^l is smallest. If is the nearest grid point above 
z and Zjjj the nearest grid point below z, then = z^ if (Zp-z) 
^ h'(Zp) < (z-Zjjj) X h'(Zj^), or z^ = z^^ if (z-zj x h'(Zm) < (Zp~z) x 
h' (Zp) . Note that h' is the inverse of the mapping metric given by 
Eq. (3.3.10). 

3.7 Generating a Smooth, Divergence-Free Initial Disturbance Field 

Kwak et al. (1975) developed a scheme for producing a divergence- 
free flow field with any desired three-dimensional energy spectrum. 
This scheme uses random numbers to select the phases of the Fourier co- 
efficients of the velocity field on a (16) uniform cubic grid. The 
Fourier coefficient vector is in the plane whose normal is the wave- 
vector, but with random phase. The orthogonality of the Fourier coef- 
ficient and the wavevector ensures that the field is divergence-free. 

We used Kwak's routine with three different sets of random numbers 

(designated 1,2,3), and two different spectrum shapes were used. In all 
• ^ 

but Case 11, we used the homogeneous Isotropic turbulence spectrum of 
Comte-Bellot and Corrsln (1971) and applied the Gaussian filter (2.5.4) 
to give a spectrum we could adequately resolve. For Case 11 we used a 
"white noise" energy spectrum (all modes excited at equal energy) as 
input to Kwak's routine. We then zeroed the Fourier coefficients of all 
modes with jkj^l >4 or |k 2 l >4. Thus we retained only 24 of the 
longest wavelength modes (in horizontal planes) in our initial distur- 
bance field in Case 11. 

Kwak's scheme generates (16) disturbance fields on a uniform cubic 
grid. From these fields we extract five planes of data and deposit them 
on the five central planes of the non-uniform anisotropic grid used in 


See Chapter 5 . 
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the mixing layer computations <> On the new grid, the disturbance field is 
discontinuous and no longer divergence- free . We first eliminated the 
discontinuity by applying a Gaussian convolution filter, as given by 
Eqs. (2.5.3) and (2.5.4) in the X 3 direction (across the mixing 
layer) . Note that this smoothing is applied to the initial conditions 
and therefore there are no difficulties of the kind discussed in Section 
2.5. The filter width used was = v'W h^ , where h 3 is the computa- 
tional grid spacing (A 5 ) . Thus the filter is non-uniform in physical 
space. We then multiplied each velocity component u^^ by a weight 

factor aj^ such that the final field has approximately equal rms amp- 
litudes for all three components . 

Finally, we obtained the divergence-free initial disturbance field 
Uj^ by two different procedures. The first, used for all but Case 11, 
was to take the curl of smoothed and weighted field. In Case 11 we 

again took the curl, but called this the disturbance vorticity and 
solved the Poisson equation for the vector potential whose curl is then 
the divergence-free velocity field which is then used in the calcula- 
tion. 
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Chapter 4 
METHOD VALIDATION 


4 «1 Convection 

In this section we shall assess the new Fourier method described In 
the last chapter by applying It to a one-dlmenslonal convection problem. 
We shall compare It to the finite-difference methods given In Appendices 
A..1 and A. 2. The solutions are all computed on grids of 33 points and 
compared with the analytical solution. 

The grid used Is defined by Eq. (3.3.9) with 


= 32/tt and Cj = (j-l)/64 


The one-dlmenslonal convection equation to be used as a test Is : 


3u , 3u „ 

IT + •= 5? ■ “ 


(4.1.1) 


and has the exact solution 


u(z,t) = f(z-ct) (4.1.2) 

which just says that any Initial waveform propagates toward Increasing 
z with 61 uniform speed c. We used a Gaussian waveform Initial condi- 
tion, for which the exact solution to Eq. (4.1.1) Is: 

u(z,t) = exp{ - [1.6651(z-ct)/6j^^2]^^ (4.1.3) 

In the problem we solved, we took ^ 1/2 ~ ^ c => -1. The half- 

wldth the width of the waveform at half Its peak value. 

We used the mapping (3.3.9) which Is appropriate for functions 
which vanish at ± Infinity. Although the appropriate Fourier method for 
this problem Is the complex exponential transform, we used an equal 
combination of sine and cosine transforms. 
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Another way of looking at this is to use the Riemann sheet perspec- 
tive of Section 3.3. A cosine expansion implies a waveform imaged sym- 
metrically about infinity onto a second Riemann sheet. Similarly, using 
a sine expansion Implies an antisymmetric reflection of the physical 
waveform onto the second Riemann sheet. The physical and image wave- 
forms will propagate towards one another with equal speed and meet at 
infinity at infinite time. If we use the combination 


6u 

1 

6u 

6u 


■5? - 

2 

•SF 

6z 

sin 

cos 


no image waveform appears on adjacent Riemann sheets , and we obtain 
periodicity of the waveform over two Riemann sheets. All waveforms are 
identical and propagate with equal speed in the same direction. 

The time advance method was the well-known fourth-order Runge-Kutta 
scheme with very small time step. The Courant number was taken to be: 


At Id 

Az . 
min 


.01 


The time step was chosen small so that the error is dominated by that of 
the spatial-differencing scheme. 

“§C 

The scheme given by Eq . (A.l) is second-order in physical space 
and slightly more accurate than Eq. (A.2) , which is second-order in 
computational space but first-order in physical space. (In the limit Aj 
^j-l> (A*2) is also second-order in physical space.) 

Figure 4.1.1 shows the grid points relative to the initial waveform 
and also shows the computed solutions at T = ct/6j^y2 °° 2.0 obtained 
using Eq . (A.l) (denoted F.D.) and the new Fourier scheme (denoted 

N.F.). The maximum error in the solution of the one-dimensional wave 
equation is .34 using Eq . (A.l), .40 using Eq . (A.2), and .0032 for the 
new Fourier scheme. Therefore, we conclude that the new Fourier method 
is vastly superior to finite-difference calculation in its handling of 
convection. 


* 


See Appendix A, Eq . (A.l). 


34 



4.2 Diffusion 


We shall now assess the new method on a diffusing problem. Using 
the same grid as for the convection problem, we solved the heat equation 


3u 3^u 

^ 


(4 .2.1) 


The fourth-order Runge-Kutta scheme was used for the time advance 
(the time was small enough that spatial differencing errors dominated) . 
The dimensionless time step was 


vAt 


.05625 


Az 


min 


We used two finite-difference schemes for the second spatial deriv- 
ative. Idle first scheme is two consecutive applications of (A.l); the 
second finite-difference scheme is given by (A.3) . The (A.3) scheme is 
a three-point scheme and is first-order; it becomes second-order as 

The initial condition used was an error function, giving the ana- 
lytical solution: 


u(z ,t) 


erf 


z/AvtJ 


We set = 25 and v = .06 , and we advanced the computation until 

tf = 16t^. 

Figure 4.2.1 shows the time history of the maximum normalized error 
defined by 


m 


(u -u.) 
c A max 

"(u '-u'l ) 

o A max 


as a function of dimensionless time. E is the maximum error in the 

m 

computed solution normalized by the maximum change from the initial 
condition. We see that the new Fourier scheme has errors six orders of 
magnitude smaller than the finite-difference method at early times and 
three orders of magnitude smaller at later times. Thus the new Fourier 
method is far superior in its treatment of diffusion. 
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4 .3 Taylor-Green Problem (2-D) 

The previous two tests were both on one-dimensional problems and 
pertained to the infinite domain Fourier method. In order to test the 
Navier-Stokes aspects of the final program, we checked the accuracy of 
the code by computing the solution of the 2-D Taylor-Green problem and 
comparing it with the analytical solution. The 2-D Taylor-Green problem 
is a stable configuration of counter-rotating vortices whose amplitude 
decays by viscous effects. The initial condition for this problem is 

u(x,y,z,0) = cos(k]^x) sin(k 2 y) 

v(x,y,z,0) =■ ~(kj^/k 2 ) sin(kj^x) cos(k 2 y) (4.3.1) 

w(x,y,z,0) = 0 

The analytical solution to this problem is 

-(kj+k2)vt 

u(x,y, 2 ,t) = u(x,y,z,0) e 

-(kj+k^)vt 

v(x,y,z,t) = v(x,y,z,0) e 

w(x,y,z,t) = 0 

In our calculation we used v = 1.36, kj^ = k 2 = ir/39.6, and At = 
6.26, which gives 

-(kJ+k^)vAt 

e = .8983775 

The fourth-order Runge-Kutta method gives the first five terms of a Tay- 
lor expansion of the exponential, giving an error of .0000011 by analy- 
sis. Exactly the same error was found in the computation. 

Since the error in the computation matches the analytical error es- 
timate , we conclude the code is functioning properly in two dimensions . 
This test showed that the pressure, convection, and diffusion are 
properly advanced in time , at least for two-dimensional flows . 
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Chapter 5 
RESULTS 


5 .1 Overview of Effects Studied 

The primary focus of this study concerns the sensitivity of the 
mixing layer to Initial conditions . We shall examine the growth of the 
momentum thickness , mean velocity profile , kinetic energy , Reynolds 
stresses, anisotropy of the Reynolds stresses, shear stress correlation 
In the center of the layer , and particle track pictures . These data 
V7111 allow us to assess the mixing layer's sensitivity to the amplitude, 
spectrum shape, and relative phases of the Initial disturbance field. 
VJe shall also examine the Influence of the subgrid scale model on the 
computed large-scale field, as well as the Influence of filtering. It 
must be remembered that, due to computational limitations, all computa- 
tions cover only the developing near field of the mixing layer . 

5 .2 Case Descriptions 

Nine cases of mixing layer transition were studied In detail . In 
Sections 2.3 and 3.7 we discussed the way In which the Initial condi- 
tions were constructed. In this section we shall give a detailed case- 
by-case description. Thus the various cases will be freshly in mind 
later In the chapter, allowing a better understanding of the results. 

The Initial conditions of the various cases are defined by their 
disturbance fields, as the initial mean field is the error function in 
all cases. Each Initial disturbance field Is characterized by the 
energy spectrum and random number set used by Kwak's routine, the method 
by which It Is adapted to the nonuniform grid, and the amplitude. The 
description for each case Is given In Table 5.1. The terms are defined 
In Table 5.2. 

In addition to the effect of Initial conditions, we also checked 
the sensitivity of our results to the subgrid scale model described in 
Section 2.6 and the Influence of the filtering the governing equations. 
The subgrid scale model is given by Eq. (2.6.1) with the eddy viscosity 
given by (2.6.2). Whether the model and/or filtering were used is 
Indicated in the last two columns of Table 5.1. Case 9 was an attempt 
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to resolve the full field without filtering, but the resolution measures 
described In Section 3.5 suggest that, by dimensionless time T = 130 
or 0/0Q = 2.4, the resolution Is marginal. The Initial condition of 
Cases 5, 9, and 10 are Identical. Comparing Case 5 with Case 10 shows 
the model Influence, while Cases 9 and 10 assess the Influence of 
filtering . 


Table 5.1 

Descriptions of Cases Run 


Case No . 

Input 

E(k) 

Amplitude 

Random 
No. Set 

Initial 
Field Type 

Model 

Filtering 

2 

C-B-C . 

High 

1 

1 

Yes 

Yes 

4 

" 

Low 

1 

1 

If 

M 

5 

t* 

Medium 

1 

1 

ti 

«• 

6 

tt 

Low 

2 

1 

” 

M 

7 


Medium 

2 

1 

M 

tt 

8 


High 

2 

1 

tl 

tt 

9 


Medium 

1 

1 

No 

No 

10 

M 

Medliam 

I 

1 

• ( 

Yes 

11 

Box 

Medium 

3 

2 

(1 

tt 


Table 5.2 

Definitions of Terms Used In Table 5.1 


Quantity 

Descriptor 

Significance 

E(k) 

C-B-C 

Comte-Bellot & Corrsln spectrum. 


Box 

White Noise spectrum. 

Amplitude 

Low 

- 3.2 X 10-6. 


Medium 

«iV2W 3.2 X 10-‘. 


High 

= 3.2 X 10-2* 

Random No . Set 

1,2,3 

Seed for random number generator . 

Initial Field Type 

1 

Curl used as Initial velocity. 


2 

Curl used as Initial vortlclty . 

Model 

Yes , No 

Subgrid scale model used? 

Filtering 

Yes, No 

Filtering used during simulation? 
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Before discussing the results , we shall comment on the length of 
the various runs. Initially, all runs were given enough computational 
time to ensure at least a fourfold thickening of the layer. Since the 
time step was varied continuously in all cases, the final times of the 
runs differed significantly. Following this initial set of runs, three 
cases ('),6,7) were run until the layer was at least eight times its 
initial thickness. We thereby explored several effects, using a few 
runs to the limit of our grid without excessive computer cost. (Each of 
the various cases required from 1.5 to 5 hours of CDC 7600 computer 
time .) 


5 .3 Mean Velocity Profiles 

The mean velocity normalized by the velocity difference across the 
layer was plotted at various dimensionless times T = AUt/6^ against 
the scaling variable Z = z/0, where 0 is the momentum thickness 


■/:[■'• -M] 


(5.3.1) 


where < > denotes a planar average. As nearly all of the profiles 
were found to be self-similar throughout the simulations, we show only 
one case here. Case 6, the case that departs most significantly from 
self-similarity . 


Figure 5.3.1 presents the mean velocity profile in terms of the 
non-dimensional variables defined above for Case 6, a low-amplitude 
initial field case which was run for a long time. It is clear that the 
mean velocity profile remains self-similar for a long time — a non- 
dimensional time of the order of 400. At the last time shown (T = 
554) , self-similarity breaks down, and this is a clear indication that 
the numerical simulation is no longer faithful to the physics of the 
flow. Ihe calculation clearly has to be stopped at this point. Mote 
that the layer was more than ten times its initial thickness at this 
time . 
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5 .4 Momentum Thickness 


In this section we shall discuss the growth of the momentum thick- 
ness. Before looking at the ntunerical simulations, we shall discuss 
some experimental results. To compare the experimental spatial growth 
rate, d6/dx, with the growth rate of the time-developing layer, 
(d6/dt)/AU, we Introduce the parameter 


U1-U2 

U1+U2 


(5.4.1) 


where is the high speed and U 2 is the low speed in the experi- 
ment. By Taylor’s hypothesis, the time rate of growth of the momentum 
thickness in the time-developing mixing layer is related to the spatial 
growth rate by : 


d6/dt 1 d6 
Au ° 2 X 3x 


(5 .4.2) 


The experimental d0/dx are converted using (5.4.2) to allow com- 
parison with the computational results. For a mixing layer to be self- 
similar, 6 must grow linearly in time or space, but the shear layer 
need not be self-similar except in the developed, far-downstream state. 
Whether there is a unique state for all shear layers is not known. 


Mansour et al. (1978) in Table 1.1 gave a list of the growth rates 
for several different experiments. The experimental data for the 
momentum thickness were fit with straight lines to give values of the 
growth rate; the resulting values of (d6/dt)/AU vary between 0.015 and 
0.022. We cannot expect a precise comparison with the experimental data, 
due to differences between a time and spatially developing layers, as 
discussed in Section 2.2. The differences are likely to be greater in 
the early or near field. In the spatially developing case, there is a 
feedback due to the pressure that is absent in the time-developing 
case. The measurements of Wlnant and Browand (1974) gave (d9/dt)/AU = 
0.035 in the near field, while farther downstream (d6/dt)/AU = 0.019 
for the experiment with laminar boundary layers at the edge of the 
splitter plate. 
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Figures 5.4.1 through 5.4.3 show ®/6 q versus T = tAU/6^ for the 
various cases . We shall now discuss the effects of Initial conditions , 
the subgrid scale model , and filtering on the development of the momen- 
t<jm thickness . 

A. Initial Conditions 

First we shall examine the effect of the different random phases of 
the Initial disturbance fields. The different phases were obtained by 
starting the random number generator with a different "seed." The small 
initial disturbance cases (4 ,6) differ only in their initial random 
phases. Figure 5.4.1 shows that the time required to attain a signifi- 
cant growth rate can vary, but Fig. 5.4.4 shows that, once a significant 
growth rate Is attained, the Initial random phases have little effect. 
A similar behavior Is observed In Cases 5 and 7, whose Initial dlstur- 
bance fields are 10 times as energetic as Cases 4 and 6, respec- 
tively, but are otherwise Identical. We again see in Fig. 5.4.2 that 
the case (7) with seed #2 grows sooner than the case (5) with seed 
#1 . This suggests that the random phases produced by seed #2 are more 
closely aligned with the phase distributions of the most amplified 
eigenfunctions. Moving to Figs. 5.4.3 and 5.4.6, we again Increase the 
initial disturbance energy by 10*. In these two high-amplitude Initial 
disturbance cases (2 and 8) , the Initial phases have a very minor role 
with seed iH , once again showing a higher early growth. 

In conclusion, we note that the Initial random phases affect early 
growth of the momentum thickness , with diminishing Influence on later 
growth. We also note that the influence is greatly diminished for high 
Initial disturbance amplitudes . 

We now examine the effect of the spectrum shape of the initial con- 
dition on the growth of the momentum thickness. Cases 10 and 11 differ 
In the spectrum of the Initial field with the energy content of Case 11 
concentrated In the amplified and most slowly decaying modes relative to 
Case 10. (Note that both Cases 10 and 11 are run without a subgrid 
model .) 

As Figs. 5.4.2 and 5.4.5 show, the Initial concentration of energy 
In the most amplified modes results In a significantly higher momentum 
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thickness growth rate. Though the present computations do not establish 
the duration of this higher growth rate, it is fair to say that a sig- 
nificantly higher growth rate due to concentrating the initial distur- 
bance energy in the most amplified modes will persist for at least a 
sixfold thickening of the layer for the case of a medium amplitude 
initial disturbance. 

Turning our attention to the effect of the initial disturbance 
amplitude, we now examine Figs. 5.4.7 and 5.4.8. Cases 4, 5, and 2 
(Fig. 5.4.7) all have the same initial spectrum shape and random phases, 
but differing initial disturbance amplitudes. We note that the small 
and medium amplitude Cases 4 and 5 display nearly the same behavior as 
each other; the low-amplitude case is slightly delayed in its growth. 
The high initial disturbance amplitude (Case 2) shows a differing trend. 
It first reaches a growth rate two to three times as large as the asymp- 
totic value of the spatially devloping mixing layer, and then plummets 
to the observed range of the asymptotic spatially developing layer. The 
behavior of Cases 6, 7, and 8 is similar to the behavior observed in 
Cases 4, 5, and 2, respectively. 

We conclude that, at low initial amplitude, different initial amp- 
litudes lead to the same momentum thickness growth rate; at high initial 
amplitudes, the growth rate may dramatically overshoot its final value. 

B . Subgrid Scale Modeling 

The effect of a subgrid model is seen by comparing Cases 5 and 
10. In Fig. 5.4.2, one sees that the model slows the growth of the 
momentum thickness only very slightly. Figure 5.4.5 emphasizes that the 
effect is indeed slight, even at later times. Thus the effect of the 
subgrid scale model is not very important in simulation of these flows; 
more evidence will be given later. We shall also see later that the 
differences that do exist are in the small scales, as expected. 

It appears that the model has little influence on the growth of the 
momentum thickness . However , the observation may be due to cancellation 
of two effects; the model destroys the small scales, making the flow 
less energetic, but it also increases the diffusion of momentum. 
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C. Filtering 


Next, we look at the effect of filtering by comparing Cases 10 and 
9. For Case 9, which has no filtering, the growth rate of its momentum 
thickness is quite erratic after T = 130 or 6/6^ = 2.5. As we shall 
see later, this case is behaving unphysically after this point. 

Thus filtering the field is very effective in keeping energy from 
piling up at the high wavenumbers. Indeed, filtering is more important 
than the subgrid scale model in this respect. 

5 .5 Turbulent Kinetic Energy at the Center of the Mixing Layer 

To examine the energy of the disturbance field, we have plotted the 

turbulent kinetic energy, ~ '^i ~ ^ '^i center 

of the mixing layer, normalized by the square of the velocity difference 

2 

across the layer (AU) , against both dimensionless time T and dimen- 
sionless momentum thickness. These curves display some important fea- 
tures that we shall now discuss . 

An interesting effect is a tendency for the (normalized) energy to 
overshoot the value expected based on the experiments . Fig. 5.5.1 shows 
the time history of the normalized kinetic energy at the center of the 
mixing layer, for Cases 4 and 6, the small-amplitude initial distur- 
bance cases. While the overall behavior of the two cases is similar, 
they begin to diverge somewhat around T = 50, but later rejoin. By T 
= 300, both are very near 0.037, the far-downstream value of the nor- 
malized kinetic energy reported by Wygnanski and Fiedler (1970). How- 
ever, both cases exhibit overshoot; Case 6 reaches a peak about 60% 
higher than the final experimental value. 

Although the calculations cannot be carried sufficiently far to 
reach the asymptotic state of the mixing layer, we believe this over- 
shoot is real, because it has also been seen in experiments. Bradshaw 
(1966) reported an experiment with a shear layer generated from a lami- 
nar boundary layer . In his results , the gradient component of the 
disturbance field < w > overshoots its far-downstream value by 100% 
(rms) and the streamwise component overshoots by 15%. His results are 
shown in Figs. 5.5.2 and 5.5.3. Thus the overshoot we find is consis- 
tent with the experimental data. 
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The computed results appear to have a longer overshoot time scale 
than the experiment; this may be due to the difference between the 
time- and space-devloping layers, or it may be a Reynolds number effect. 
Although not shown, the kinetic energy of Case 6 begins to decline at T 
= 450 and continues to descend toward the far-downstream value until 
the computation is stopped at T = 554. However, the resolution is 
marginal from T = 450 on. In Fig. 5.5.4, the normalized kinetic energy 
of Cases 4 and 6 is plotted versus the dimensionless momentum thickness 
0 = Q/Oq* these coordinates the two cases are remarkably similar, 
and this indicates that the increase in kinetic energy is more related 
to the thickness of the mixing layer than to the time. Hence, we shall 
plot (q^/2)/(AU)^ vs. 0/0Q in the remainder of this section. 

Let us now look at the overshoot in the medium amplitude cases. 
Figure 5.5.5 shows the normalized turbulent kinetic energy as a function 
of dimensionless momentum thickness for the medium- amplitude disturbance 
cases. The initial disturbance kinetic energy is slightly more than two 
orders of magnitude smaller than the far-downstream turbulent kinetic 
energy. In the two cases (5 and 7) run furthest in time, we see an 
overshoot followed by a gradual decline. The overshoots are approxi- 
mately the same size as in the low Initial disturbance cases. Case 9 
differs most from the rest of this group . This is the case with no 
filtering or subgrid model. After reaching a peak amplitude roughly 45% 
higher than the far-downstream experimental value , the kinetic energy 
begins to decline. Unfortunately the numerical resolution of this case 
is questionable after T = 130 (0= 2.4). 

A striking observation is that the overshoot behavior at high ini- 
tial amplitudes differs from that expected from experiments . The two 
high-amplitude cases (2 and 8) are shown in Fig. 5.5.6; these cases 
might represent the mixing layer produced from turbulent boundary lay- 
ers . In these cases the initial energy is roughly 10% lower than the 
far-downstream experimental value. Both of these cases show a quick 
rise of the kinetic energy to roughly 60% more than the far-downstream 
value of the laboratory layer; at the time of the peak in the energy, 
the mixing layer has thickened by a factor of 2.5. Following this, the 
turbulent kinetic energy for Case 2 decreases , but it is still about 35% 
above the fully developed turbulent mixing layer laboratory value at the 
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end of the run. Case 8 shows no sign of a decline In the turbulent 
kinetic energy at the time the run was stopped . 

The most significant aspect of these two hlgh-amplltude runs Is 
that the energy and growth rates significantly overshoot the expected 
fully developed values . However , the experimental mixing layer does not 
exhibit this behavior when the splltter»plate boundary layers are turbu- 
lent. Tlie most probable reason for this discrepancy Is that the bound- 
ary layer turbulence has a longer characteristic streamwlse wavelength 
and much higher obliqueness (Kim, 1981) than the mixing layer can amp- 
lify or sustain, while the Initial disturbance used In the simulation 
contains more of the highly amplified wavelength components . We must 
also consider the possibility that the time-developing mixing layer is 
more energetic than the spatially developing layer . 

The effect of the subgrid model on the centerline kinetic energy is 
seen by comparing Cases 5 and 10 ,whlch are Identical In all aspects , 
except tliat no model is used in Case 10 in Fig. 5.5.5. The model re- 
sults in roughly a 10% reduction in the energy, which is smaller than 
the effect produced by varying the initial conditions (compare Cases 5 
and 10 with Cases 7 and 11) . 

5 .6 Profiles of the Turbulent Kinetic Energy 

In the previous section the turbulent energy at the center of the 
layer was presented. That Information is sufficient to gain an overall 
perspective of the turbulence Intensity. In this section, we shall 
check for similarity in the turbulent energy profile. 

In Fig. 5.6.1, the turbulent kinetic energy, normalized by the 
value at the center of the layer, is plotted versus distance across the 
layer normalized by the momentum thicknesses for Case 4 , a low Initial 
amplitude case. Case 6 Is quite similar and Is not shown. In the time 
Interval covered in the figure, initial momentum thicknesses grew by 
between 3 .5 to 4 .0 . 

In the small disturbance cases, the Initial profile was much too 
broad. At T = 83 we find a narrow profile characteristic of the 
eigenfunctions of the linearized equations. The profile then begins to 
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broaden until, at T = 246, the energy profile is somewhat broader 
than the "fully developed turbulence" case of Spencer and Jones (1971) . 
At the latest time, the profile is close to that of Spencer and Jones. 
To see whether this profile persists, we ran Case 6 out much further in 
time. The profiles remain close to the experimental one; at the final 
time the layer is more than ten times as thick as it was initially and 
the numerical accuracy is deteriorating, due to the coarseness of the 
grid at the outer edges of the shear layer. 

In Fig. 5.6.2, the normalized turbulent energy profile is given for 
a typical medium-amplitude case. Case 11; the others are similar. We 
again note that the initial condition profile was much too broad. At 
approximately T = 66 , these cases exhibit the "narrow" profile char- 
acteristic of the eigenfunctions of the linearized equations. At later 
times all cases show profiles characteristic of fully developed turbu- 
lence without the oscillatory behavior of the profile width found in the 
small-amplitude cases . 

In Fig. 5.6.3, we give the kinetic energy profiles of a high 
initial- amplitude case (Case 2); Case 8 is quite similar. In these 
cases the profile exhibits a self-similar shape by T = 32.8. While the 
shape reaches self-similarity very quickly, we must remember that the 
maximum intensity overshoots the final value due to the long wave char- 
acter of the initial condition. 

In summary, we note that the profile of the turbulent kinetic en- 
ergy profile always progresses from the broad initial profile to the 
self-similar experimental profile. The approach is oscillatory when the 
initial energy is low and mono tonic if the initial energy is high. 

5 .7 Reynolds Stresses 

In the previous sections, we have learned something of the behavior 
of the mean velocity profile, the momentim thickness, and the turbulent 
kinetic energy. With this in mind, we shall now look at profiles of the 
Reynolds stress tensor. The discussion will proceed from the small- 
amplitude to large-amplitude initial disturbances with the effects of 
modeling and filtering included in the medlumramplitude initial distur- 
bance cases . 
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In Figs. 5.7.1 through 5.7.7, the profiles of the Reynolds stresses 

for a small- amplitude case (6) are presented. Because the absolute 

values of the quantities increase rapidly, the results are presented as 
~ ~ 2 

< >/q (z=0) vs. z/6. The first figure shows the initial 

conditions . The maximum value of each normal stress is nearly the 
same. The double peaks in the normal stresses < u > and < v > are 
due to the way the initial condition are constructed. 

At T = 83 (Fig. 5.7.2), the stresses are those characteristic of 

the eigenfunctions of the linearized equations. By T = 165, nonlinear 

^2 

effects are becoming important and the dominance of < u > is dimin- 
ishing; see Fig. 5.7.3. At T == 246 (cf. Figs. 5.7.4), we see that the 
gradient component of the normal stress , < W^ > , has become strongly 

dominant . 

A physical explanation of this can be offered with the aid of the 

visualizations that are presented in Section 5.11. At T = 246, the 

shear layer has rolled up into vortical structures. If these were 

straight two-dimensional vortices, as Brown and Roshko (1972) suggested 

~2 ~2 

they might be, < u > would be double peaked, < w > would be single 

~2 ~2 

peaked and larger than < u > in the center, and < v > would be 

zero. This picture explains much of what is seen in Fig. 5.7.4. The 

differences are due to the vortices not being straight, which introduces 

~2 

streamwise vorticity, causing < v > to be non-zero, and reduces 

< u >. Bradshaw (1966) found the same effect in the laboratory. In 
fact, his data show a stronger effect; in some cases the gradient stress 
was three times the streamwise stress. 

At T = 328, Fig. 5.7.5 shows the onset of significant asymmetry. 
Tills is not surprising as there are only two or three large eddies in 
the computational domain. The strong dominance of the gradient compo- 
nent of the normal stress is greatly diminished at T = 328. Figures 
5, .7. 6-7 show later profiles for Case 6. Although the streamwise compo- 
nent becomes the dominant normal stress at T = 554 , the results are 
marginal at this time due to Insufficient resolution at the edges of the 


The initial field is defined by x v^, v^ varies rapidly 

with z near the edges of the layer. 
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shear layer. The double peak of the streamwlse component is probably 
due to the existence of a strong vortex with some curvature in the plane 
of homogeneity (see Section 5.11). 

Finally, we note that < w^ > is dominant at the outer edges of the 
layer at late times but not in the initial condition. This is consis- 
tent with experimental observations and, as Phillips (1954) has shown 
~2 ~2 ~2 

analytically, <w > = <u > + <v >in the irrotational portions of 
the flow. This remark applies to other cases as well. The other small- 
amplitude case (4) is similar to Case 6, but most of the effects ob- 
served above are weaker. 

Let us now turn our attention to the medium- amplitude cases (5, 7, 
9, 10, 11). The initial conditions are shown in Fig. 5.7.8-10. It is 
apparent that the method used to generate Case 11 results in much 
smoother stress profiles and probably should be preferred. 

Figures 5.7.11-14 show four of the profiles at T = 66. At this 
time the profiles have the characteristic shape corresponding to the 
eigenfunctions of the linearized equations. All cases are now very 
similar, in contrast to the initial conditions. The S 3 raimetry is appar- 
ent , as is the dominance of the streamwlse normal stress . Cases 5 and 
10 (not shown) are nearly Indistinguishable, showing that the subgrid 
model has little Influence on the stresses; this accords with what we 

found earlier about the need for the model. The "full simulation" (Case 

^^2 ^2 

9) yields rounder profiles for < v > and < w > than the filtered 
calculations . 

In Figs. 5.7.15-18, we see immense differences between cases devel- 
oping as we move into the strongly nonlinear domain at T - 131. Cases 
5 and 10 (not shown) , which differ only by the absence of a subgrid 
model in Case 10, again display very little difference. However, 
comparing Cases 5 and 10 with Case 9, we note a large influence of 
filtering. At T = 130, Cases 5 and 10 are only roughly half as 
energetic as Case 9 at T = 130. At T - 180-190, Cases 5 and 10 are 
as energetic as Case 9 at T “ 130. If Fig. 5.7.17 is compared with 
5.7.19, it is found that they are quite similar. This suggests that the 
influence of filtering is to delay the time required to reach the 
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strongly nonlinear domain. Case 9 is generally closer to the physics of 
the laboratory flow, as might be; expected. 

The difference between Cases 5 and 7 at T - 131 is amazing, since 
the two cases differ only by the set of random numbers used to generate 
the initial conditions. Figs. 5.7.19-22, give the stress profiles at T 

- 197. Cases 5 and 10 are again nearly Identical. Cases 5 and 7 are 
much more similar than at the previous time but, as we shall see below, 
this is probably coincidental. Case 9 is significantly different from 
Cases 5 and 10, but Case 9 has marginal resolution at T = 197, so cau- 
tion is needed. The breakdown of the resolution is hinted at by the 
jaggedness of the curve; better evidence of breakdown will be given 
later, cf . Sections 10 and 11. Case 11 is now quite different from the 
others. Once again, this demonstrates the ability of nonlinear pro- 
cesses to cause a large divergence of solutions which were close to one 
another (compare Cases 10 and 11 at T =* 131 with the same cases at T 

- 197). Figure 5.7.23 shows that Case 5 has nearly reached the asymp- 
totic state at T “ 295; no further significant changes were observed 
at later times. On the other hand. Figs. 5.7.24-25 show that Case 7 
takes a long time to reach the far-downstream state. In Section 5.8 we 
shall present stronger evidence of this. 

Next, we shall compare the medium Initial amplitude cases with the 
corresponding low initial amplitude cases; Cases 4 and 5 and Cases 6 and 
7 have the same initial fields except for the magnitude of the distur- 
bance field. Comparing the latter pair of cases, we find very similar 
behavior, but the medium Initial case (7) develops on a faster time 
scale than the low initial amplitude case (6) . The same is true of 
Cases 4 and 5 . 

Let us now turn to the high initial amplitude cases (2 and 8) . The 
initial profiles are the same as those for Cases 4 and 6, respectively. 
At T - 33, Figs. 5.7.26-27 show that in Case 2 the streamwise nor- 
mal stress is more dominant than in Case 8. both of these cases show 
stronger dominance of the streamwise component than that observed 
experimentally with turbulent boundary layer(s) on the splitter plate; 
they also differ from the experiment in other respects, as we have seen 
earlier and shall see again. Figures 5.7.28-30 show stress profiles of 
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roughly the experientally observed relative magnitudes , indicating that 
these cases develop more quickly than the low amplitude ones . Figure 
3.7.29 shows that the gradient normal stress becomes dominant for Case 8 
at T = 66; this is not observed in experiments having turbulent bound- 
ary layers. Fig. 5.7.30 shows that Case 2 has considerable asymmetry in 
the spanwlse normal stress; this reflects the small sample size. All of 
the differences between our computations and experimental results for 
the high initial amplitude cases (except the asymmetry) show the compu- 
ted results to be in the direction of the experimental results for 
laminar boundary layers . 

We thus conclude once again that the experimentally observed near- 
field differences between mixing layers generated from laminar and tur- 
bulent splitter plate boundary layers are due more to the characteristic 
wavelengths of boundary layer turbulence than they are to the distur- 
bance amplitude. Additional conclusions about the Reynolds stress evo- 
lution will be drawn at the end of the next section. 


5 .8 The Reynolds Stress Anisotropy Itensor 


In this section we shall discuss the Reynolds stress anisotropy 

/V o/ rv 

tensor. We define < u.u . > as the average of u^u. over a hori- 

i 3 i J 


zontal plane . 


If there were enough points in each horizontal plane , 


< > would be equivalent to the ensemble average of u^u^.i.e., 
the Reynolds stress tensor. However, as there are only 256 points in 
each horizontal plane and the velocities are not statistically indepen- 
dent, this identification needs to be made cautiously. The anisotropic 
component of this "Reynolds stress" is : 


'IJ 


< > 


3 *^ij 


(5.8.1) 


In the time developing mixing layer, b^j is a function of z 
and t. As the Reynolds stresses are functions of z and there are 
asymmetries , we shall consider an energy-weighted average of the 
anisotropy tensor (5.8.1) to characterize the turbulence 
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lliis will limit the number of figures required to present the results 
and carries most of the important information. 

We shall begin examining the energy-weighted anisotropy tensor by 
looking at the small-amplitude cases. Wext , the medium-amplitude cases 
will be considered, and here the effects of modeling and filtering are 
examined, as well as initial condition influence. Finally, the high 
initial disturbance amplitude cases are considered. In Figs. 5.8.1 

through 5.8.9 the for the nine cases studied are plotted versus 

time . 


The for the small-amplitude cases (4 and 6) shown in Figs. 
5. 8. 1-2 are similar, but not identical. The streamwise normal stress 
dominates at early times , while the gradient component dominates at 
Intermediate times. The shear (B13) and span (B22) ts have similar 
character, and both are primarly negative. For the shear component, 
this is as expected. The span component is the smallest of the normal 
stresses , which is different from most observtions of the asymptotic 
state. Even at T = 550, at which time the layer is 10 times its ini- 
tial thickness , the anisotropy tensor of Case 6 is far from the asymp- 
totic one. 

Next, let's look at the medium-amplitude cases. The Bj^j for the 
medium-amplitude cases (5,7,9,10,11), are shown in Figs. 5.8.3 through 
5.8.7. Comparison of Figs. 5.8.1 and 5.8.3, which represent cases with 
the same initial field except for the amplitude, are quite similar, 
except that in the case with higher initial amplitude (5), the gradient 
component of the normal stress never becomes dominant. The case shown in 
Fig. 5.8.4, which is a higher amplitude version of the case shown in 
Fig. 5.8.2, does display dominance of the gradient component at later 
times , but the dominance is not as strong as in the low initial distur- 
bance level case. Thus is appears that increasing the amplitude of the 
initial disturbance decreases the dominance of the gradient component of 
the normal stress; this will be further verified when we consider the 
high amplitude cases below. 
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Figures 5.8.3 and 5.8.6 represent two cases which are identical 
except that the latter one does not have a subgrid scale model. As can 
be seen, the effect of the model on is quite small; this is consis- 
tent with what we have observed about the effect of the model earlier. 
The case (9) shown in Fig. 5.8.5, which is identical to Case 5 just 
considered except for lack of filter and model, is similar to the other 
two; the most significant difference is that the spanwise normal stress 
plays a more important role at earlier times in the unfiltered case. 
This behavior is closer to what has been observed in the laboratory than 
the other two cases. The shear stress is also smaller and closer to the 
experimental data in this case. It seems that a significant portion of 
the spanwise stress resides in the small scales; further evidence of 
this will be given later. However, the small scales seem to destroy the 
shear stress. Unfortunately, in this case the numerical method became 
unreliable at a much earlier time, and we could not follow the devel- 
opment as long as in the other two cases with the present code. The 
solution to this difficulty lies in constructing a code with a greater 
number of grid points in the horizontal directions, and we recommend 
that this be done. The case (11) shown in Fig. 5.8.7, with initial con- 
ditions which are relatively deficient in small scales, produces greater 
dominance of the streamwlse nomal stress at early times and of the gra- 
dient stress at later times. Also, the high shear stress is maintained 
longer; this is a reflection of the more rapid growth of the layer in 
this case. 

Finally, we examine the high- amplitude cases. Figures 5.8.8 and 
5.8.9 show the for the high amplitude cases (2 and 8). The changes 
noted in going from the lowest amplitude to the medium amplitude are 
even more obvious in these cases. Relative to the case shown in Fig. 
5.8.3, which has the same initial condition except for amplitude. Fig. 
5.8.8 shows less dominance of the streamwlse normal stress at early 
times and of the gradient stress at later times; all of the stresses 
show less variation with time, and the asymptotic state seems to be 
reached more quickly. The same effects are seen in comparing Figs. 
5.8.4 and 5.8.9, but the change is perhaps less dramatic. 
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In summary, in some cases the components of Reynolds stress ani- 
sotropy tensor appear to reach asymptotic values after the layer has 
thickened only by a factor of four. In other cases no asymptotic value 
is in evidence after a thickening by as much as a factor of 10. This 
illustrates the enormous sensitivity of the Reynolds stress anisotropy 
tensor to the initial disturbance field. This sensitivity is more 
strongly coupled to the phase relationships of the initial conditions 
than to their amplitude or spectrum shape. The fact that we have only a 
couple large eddies in our computation increases this sensitivity; it 
also shows that large eddies can vary considerably, though variations 
may be small in a given laboratory experiment, due to phase locking. 


5 .9 Correlation Coefficient at the Center of the Mixing Layer 

In this section we shall examine the correlation coefficient of the 
Reynolds' shear stress. It can be defined at the center of the layer as 


c(uw) = < u(x,y,o)w(x 


,y*o>l j< Qu(x,y,o)^J > < j^(w(x,y,o))^ J >| 


Tlie correlation coefficient can range between -1.0 and 1.0. In most 
shear flows studied in the laboratory it has a value of -0.45 ± 0.05. 
However, Shlranl et al. (1981) found values as large as -0.75 in a sim- 
ulation of a homogeneous shear flow. 

The correlation coefficients for the seed #1 cases are shown in 
Fig. 5.9..1, and those for seed #2 are shown in Fig. 5.9.2. The collapse 
of the data is remarkable. llie values also seem realistic for the 
mixing layer, particularly in the near field. 

Next, we shall look at the effect of changing the seed while main- 
taining constant initial amplitude. The correlation coefficients for 
the three amplitudes used are shown in Figs. 5. 9. 3-5. For the low and 
medium amplitude cases, the collapse of the results is considerably 
poorer tlian is the case when the seed is held constant and the amplitude 
varied. The exception is the high amplitude set of cases, for which an 
excellent collapse of the data is found. 

Thus we find that the phases in the initial conditions have a 
stronger effect on the correlation coefficient than do the amplitudes in 
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the low to medium initial-amplitude cases. For high initial amplitude 
cases, the initial phases are not very Important. 

5.10. Resolution and Statistical Validity 

In Figs. 5.10.1 through 5.10.5, the accuracy parameter n^, the 
alias-free fraction and sample number Tig defined in Section 3.5 
are plotted. 

Comparison of the sample number of case 5 (Fig. 5.10.2) with that 
of case 10 (Fig. 5.10.3) shows that the subgrid scale model extracts 
most of its energy at the higher wavenvimbers , as it should. 

Recalling that cases 5, 9, and 10 all have the same initial condi- 
tion, we note that case 9 (no filtering or subgrid model), shown in Fig. 
5.10.4, shows a strong Increase in the sample number, and resolution is 
questionable after, say, T = 140. On the other hand, the Case 10 (no 
model) results are nearly the same as those of Case 5 . Filtering thus 
has a much larger Influence than does the model. 

In contrast to general decline of rig in the small and medium- 
amplitude cases. Figs. 5.10.4-5 show that in the highramplitude initial 
disturbance cases, rig hardly changes at all. In general, the effects 
of initial amplitude on these parameters are not large. 

The figures show that, with the exception of Case 9, there is good 
spatial resolution in the central horizontal plane. Again, except for 
Case 9, two to four "characteristic" large eddies are captured in the 
computations. This is a very small sample, so extreme care is needed in 
comparing the results of these computations with those of experiments 
which are averages over many large eddies . 

5.11. Particle-Tracking Visualizations 

To complement the quantitative information given in the previous 
sections flow visualization pictures are presented in this section. The 
flow visualization pictures were obtained by tracking the intersections 
of the grid shown in Fig. 5.12.1. Plan views are views along the z 
axis, with the mean flow in the horizontal direction. One can consider 
the grid lines as dye or smoke lines put into the flow. For Cases 5 
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through 1.1, this grid (5.12.1) was placed in the horizontal midplane of 
the mixing layer at the initial time. The intersections move with local 
fluid velocity; the numerical method for doing this is described in 
Section 3.6. For the small (4 and 6) and medium (5,7,9,10, and 11) 
amplitude initial disturbance cases, the vertical lines in Fig. 5.11.1 
are essentially vortex lines because the vorticity of the disturbance 
field is weak compared to the mean flow vorticity in these cases. For 
the high-amplitude disturbance cases , the vertical lines are only a 
crude approximation to vortex lines , as the initial disturbance 
vorticity is no longer negligible compared to the vorticity of the mean 
field. 

In order to help the reader understand the visualizations, we shall 
now present some basic background on shear-layer instability. The pri- 
mary instability of a time-developing mixing layer is the Kelvin- 
Helmholtz instability. The most amplified eigenfunction according to 
linear stability theory is spanwlse uniform and, if acting alone, would 
ultimately result in a vortex array which is uniformly spaced in the 
streamwise direction (see Betchov and Criminale, 1967). 

In the computations presented here, broad spectrum, random phase, 
finite-amplitude, 3-D initial disturbance fields are introduced. This 
random disturbance can be expressed as a linear combination of the 
eigenfunctions of the linearized equations and the eigenfunctions which 
grow most rapidly in time will become dominant if the initial distur- 
bance amplitude is sufficiently small. In the low amplitude initial 
disturbance cases, a number of the eigenfunctions of the linearized 
equations grow rapidly so that the flow is not dominated by the single 
most amplified 2-D eigenfunction, l.e., the 3-D eigenfunctions are 
important . 

The time development of the marker lines for Case 6 is shown in 
plan view in Figs. 5.11.2 through 5.11.7. In Fig. 5.11.2, we see the 
roll-up of the vorticity layer into two more or less coherent spanwlse 
vortex structures (similar to what the primary Kelvin-Helmholtz insta- 
bility would produce) . A mechanism of secondary instability is also 
suggested by this figure. The secondary instability is caused by the 
straining field created by the primary vortex structures. To see how 
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this occurs , 
shown below. 


let us consider the spanwlse view of the mixing layer , 


C 



The nearly two 2-D vortex structures (A + B) produce straining field or 
stagnation line flow midway between them as indicated by the arrows . 
Vortex A will entrain Irrotational fluid from quadrant A' , while B 
entrains from quadrant B' . Also note that, if the vortex structures 
have spanwlse variation in strength or position, the stagnation line 
will not be straight. A plan view of plane C-C is shown below. 



The wavy solid line represents the stagnation line, and the dashed 
line might represent a vortex line in the "braid" between vortices A and 
B. Due to the wavlness of the stagnation line, the vortex line will be 
pulled in the directions shown by the arrows . This will turn the vortex 
line, aligning it with the flow. Eventually, most of the vorticity in 
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the secondary vortices will be aligned with the stretching direction of 
the straining field. We refer to this turning and alignment process as 
"secondary instability." 

The secondary instability produces streamwise and gradient 
components of the vorticity. This completed state has been observed at 
Ciiltech by Konrad (1976) and Roshko (1980). The secondary instability 
is a very efficient mechanism for entrainment of Irrotatlonal fluid from 
the outer part of the flow into the primary vortex structure. Roshko 
(1980) shows that the secondary instability Increases the mixing in the 
shear layer and that d6/dt Increases by 10%. A similar mechanism of 
entrainment was proposed by Corcos (1981) . 

With this mechanism for formation of the secondary vortices in 
mind, let's now discuss the other cases. The discussion proceeds from 
the small to the large amplitude cases . 

In Case 6, the primary and secondary vortex structures are clearly 
delineated and are well formed by T = 200 (Fig. 5.11.3). The primary 
vortices each have slightly less than half the total circulation. The 
part of the secondary vortex marked a-b in Fig. 5.11.3 contains three 
vortex lines suggesting that it has a circulation of approximately 30% 
of the circulation of the primary vortices; this is larger than the 
value suggested by Roshko (1980) . There is also a weaker secondary 
vortex forming in the center of Figs. 5.11.2-3, but it is harder to 
see. In Figs. 5.11.4-5, the two primary vortex structures are moving 
closer together and beginning to pair, and Figs. 5.11.6-7 show a span- 
wise view of this process. It should be noted that, while the primary 
vortices are not straight, they are nonetheless pairing in an essen- 
tially uniform fashion along the length of the vortex structures. 

Interesting effects are seen near the centers of the two primary 
vortices. The one on the right hand side of these figures is more 
curved in the early stages of the rollup. Tlie streamwise component of 
the vorticity resulting from the curvature apparently enhances the 
entrainment of irrotatlonal fluid into the center of the vortex. This 
does not happen in the straighter vortex on the left. When the vortices 
begin to pair, the fatter section of the right hand vortex is heavily 
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strained and Is torn apart. This appears to be another mechanism for 
the creation of streamwlse vortlclty. 

Case 7, which Is Identical to Case 6 except that the Initial 
disturbance energy Is two orders of magnitude larger, behaves In a 
qualitatively similar way to Case 6. Fig. 5.11.8 Is quite similar to 
Fig. 5.11.4 and we see that Increasing the amplitude has caused the 
development to take place faster but has not produced any Important 
structural changes . 

Case 8, which has an Initial disturbance energy two orders of mag- 
nitude greater than Case 7, behaves In a more chaotic manner. Recall 
that In this case the vertical marker lines of Fig. 5.11.1 are not vor- 
tex lines. Fig. 5.11.9 shows that there Is a tendency for the markers 
to agglomerate In this case but the pattern Is different. The develop- 
ment Is also more rapid. 

Cases 5 and 10 have Identical Initial conditions; the only differ- 
ent Is that Case 10 has no model. They develop in a very similar manner 
and only Case 10 will be presented here. These cases differ from Case 7 
only in the set of random phases in the Initial condition. At T == 76 
(Fig. 5.11.10), it is immediately apparent that Case 10 does not roll up 
into a pair of well defined vortices. The vortices are inclined and 
branch in a complicated manner not seen in Case 6; local (rather than 
uniform) vortex pairing is taking place. The distinction between local 
and uniform pairing is probably the key to reconciling the Bradshaw and 
Browand-Roshko models of the development of the shear layer. In 
Chandrsuda et al’s (1978) visualization, there is local (In Bradshaw's 
terminology, helical) pairing, while in Winant and Browand's (1974) 
visualization there Is essentially uniform pairing along the entire span 
of the vortices. Uniform pairing may persist indefinitely, as Browand's 
experiments suggest, if three dimensional perturbations can be excluded 
In the early stages of development; the strong two dimensional vortices 
will not be modified greatly by weak three dimensional perturbations 
after rollup has been completed. On the other hand, if the near field 
has local or helical pairing , it will probably break down into the 
chaotic state normally associated with turbulence. The development of 
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Case 10 is shovm in Fig. 5.11.10 through 5.11.15. The local or helical 
pairing is immediately apparent . 

When filtering is eliminated, as in Case 9, the initial rollup is 
not very different (compare Figs. 5.11.10 and 5.11.6). However, the 
later development is much more Irregular as can be seen by comparing 
Figs. 5.11.14 and 5.11.17 and Figs. 5.11.15 and 5.11.18. It is diffi- 
cult to be certain of the reasons for this but it appears that the dif- 
ference may be due to the fact that filtering tends to remove (or, at 
least , diffuse) kinks in the vortices . Since highly kinked vortices are 
very energetic, they are capable of producing the chaotic pattern seen 
in Figs. 5.11.17 and 5.11.18. Also recall that the numerical methods 
lose their accuracy after about T = 140 in this case. 

In Case 11, the initial condition is richer in the large scales 
and leaner in the small scales, relative to the other cases. Neverthe- 
less, the development of this flow is not greatly different from those 
presented above. At the relatively early tljme shown in Fig. 5.11.19 
there are three vortices. The first two are close together at the upper 
left while the other two are close at the center right of the figure. 
It is not surprising that this arrangement leads to local pairing of the 
center vortex with the left one at the top and with the right one near 
the center. Fig. 5.11.20 clearly shows this to be the case; also note 
that the left vortex has developed a considerable amount of curvature at 
the time corresponding to this figure. Tb understand what happens next 
it is Important to remember that the boundary conditions are periodic 
and the flow actually being computed can be constructed by repeating the 
figure to the left and right (and above and below) itself to form an 
infinite array. The center of the left vortex now begins to pair with 
the center of the right one. This causes the left vortex to become 
still more curved and the resulting streamwlse vorticity pushes the part 
of the right vortex closest to it downward. This produces a downward 
kink in the right vortex and the marker points seem to be drawn together 
in the spanwlse direction when viewed from above. The result is shown 
in Fig. 5.11.21. The spanwlse view of this flow is shown in Figs. 
5.11.22 and 5.11.23. 
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The figures and arguments presented in this section provide insight 
into the nature of the processes that occur in mixing layers . The ini- 
tial layer always seems to rollup into a configuration in which the 
spanwlse vorticlty is more concentrated than it was Initially. The pre- 
cise structure of the rolled up layer seems to depend on the initial 
disturbance field in a way that would be hard to predict simply knowing 
the spectrum. It is possible to produce layers with straight two- 
dimensional vortices, ones with curved vortices, and ones with vortices 
arranged in a 'fishnet* pattern. The subsequent dynamics of the layer 
depends very much on the configuration produced by the initial rollup. 
Pairing seems to play an Important role in the later development , but it 
may take the form of either simple uniform pairing of two vortices or 
local pairing in which different parts of the same vortex pair with 
parts of different neighboring vortices. Local pairing leads to much 
more three dimensionality than does simple uniform pairing, but it looks 
as if there will always be large regions of concentrated vorticity in 
the mixing layer and that these will grow by agglomeration as the layer 
develops . 

Finally, we should note that all of these results are for the time- 
developing mixing layer. In the spatially-developing mixing layer 
studied in the laboratory there are Important feedback effects. Pres- 
sure fluctuations created in the downstream part of the flow can influ- 
ence the upstream development strongly. It is quite possible that this 
feedback can cause the layer to 'lock on' to a particular structure 
which is then maintained for a long time. 
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Chapter 6 


CONCLUSIONS 

This study of the mixing layer has focused on the physics of tran- 
sition and early turbulence in a time-developing mixing layer. The 
effort was concentrated on the effect of initial disturbance field on 
the later character of the turbulence. 

Concern about the possible influence of image flows (which are 
implied by the boundary conditions) on the computations led to the de- 
velopment of a new infinite domain orthogonal function expansion. Use 
of the discrete form of this new method in the computations eliminates 
the influence of image flows by keeping them an infinite distance from 
the flow of interest . 

This new and very accurate numerical differencing and integrating 
scheme for infinite domains was presented. It is based on the use of 
Fourier expansions and takes advantage of the computational efficiency 
of the fast Fourier transform. The new method is applicable to more 
general boundary conditions than the standard Fourier method, due to the 
use of mapping functions. (The simplest boundary conditions to imple- 
ment are periodicity, or zero, or zero-derivative conditions, or combi- 
nations thereof.) However, the allowed mapping functions are restricted 
for reasons of efficiency and accuracy. For more detail, see Chapter 3. 

Two particular mapping schemes, both for doubly infinite domains, 
were implemented. One was chosen to handle jet-type flows, while the 
other was designed for the mixing layer. Both schemes were applied to 
linear test equations having known analytical solutions . The new scheme 
was shown to have errors as much as six orders of magnitude smaller than 
common finite-difference schemes for equal numbers of mesh points . 

Using the new inf ini te- domain scheme, a 3-D, time-dependent, large- 
eddy simulation study of transition and early turbulence in a time- 
developing mixing layer was undertaken. The primary focus of this study 
concerned the effect of the initial disturbance field on turbulence 
development. Effects due to filtering and modeling were also examined. 
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To sort out the effects of the initial disturbance field, the same 
laminar, mean-velocity profile was used as the initial mean field in 
all cases. To this mean velocity field, an initial divergence-free 
disturbance field was added. We used nine cases Involving seven differ- 
ent initial disturbance fields . These seven cases allowed us to examine 
the influence of the disturbance amplitude, spectrum shape, and random 
phase sets on the resulting early turbulence. 

The computations provided the mean velocity profile, the momentum 
thickness, the turbulent kinetic energy, the Reynolds stress tensor, the 
Reynolds stress anisotropy tensor, and particle tracking pictures. 
Examination of these results provided new understanding of the mixing 
layer. Key results of this work are summarized below: 

• Self-similarity in the mean velocity profile develops very 
quickly; the self-similar profile is independent of initial 
conditions . 

• The momentum- thickness growth rate is strongly influenced by 
the initial disturbance-spectrum shape. 

• Interesting oscillatory behavior occurs in the width of the 
kinetic energy profile for the small- amplitude initial dis- 
turbances . This oscillatory behavior is not present if the 
initial disturbance is large. 

• The anisotropy tensor is a very sensitive measure of self- 
similarity. Even changing the random phase distribution in 
the initial disturbance field produces enormous differences 
in the evolution of the anisotropy tensor. 

Probably the most significant aspect of the study was revealed in 
the par tide- track pictures. Large coherent structures readily appeared 
(some similar to those of Winant and Browand (1974) and others like 
Chandrsuda et al. (1978). More Important, the mechanism for producing 
the secondary vortices was identified. These vortices develop as a 
result of spanwlse variations in the strength or position of the pri- 
mary vortex structures, which give rise to spanwlse variations in the 
straining field stagnation line. This causes the formation of pairs of 
counter-rotating secondary vortices aligned with the straining field. 
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The statistical and structural characteristics of the mixing layer 
are very sensitive to the phases of the initial disturbance. This may 
explain the differences that luive been observed among different exper- 
imental setups, each of which produces a given type of large eddy struc- 
ture. Ihus , large eddy variation may be small in a given experiment, 
but significant variations may occur from experiment to experiment. 

IWo cases were run to examine the effects of filtering and subgrid 
turbulence modeling. We found that filtering delays the onset of non- 
linear effects and gives us less than the total picture. However, it 
considereibly extends the length of time over which the computation is 
meaningful. The subgrld-scale model was shown to have very little 
influence on the calculation of the early stages of transition. 
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Appendix A 

FINITE DIFFERENCE METHODS 


For non-unlformly spaced data, we may define 


■ *}H ■ 


and use the three-point finite-difference approximation 


6u 


= auj_^ + bu. + 


(A.l) 


where 




2 > 

A + A /A 

j y j-1 


b = - a - c 


From a Taylor series, expansion (A.l) can be shown to be second-order 
accurate in Aj . On the other hand , the two-point formula 
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with 


a = - b 




is first-order accurate and becomes second-order accurate as ^ • 
For the second derivative, the approximation; 


6 u 


fix 
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with 




(A.3) 


a = 2 /[Aj_j^(Aj + A._i)] 

b = - 2/AjAj_j^ 

c « 2/Aj(Aj + 

\ 

is first-order accurate and becomes second-order accurate as ^j-i "*■ 
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Appendix B 


Given 


TIME-INTEGRATION SCHEMES 


dy/dt * f(y,t), the second-order Adams-Bashforth scheme is 




y“ + AtCl.Sf"^ - 0.5f"^^) 


(B.l) 


The fourth-order Runge-Kutta scheme is 

(B.2) 

where 

*n+-l/2 n , At n 

y = y + ^ f(y ,t ) 

y = y + y- f(y ,t ) 

- y" + it 
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The sedond-order Runge-Kutta method is 
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Appendix C 
SUBGRID MODELING 


C.l A New ttodel 

Recall from Section 2.5 that subgrid scale field Is the difference 
between the velocity and filtered velocity, denoted u^. While we 
refer to u| as "subgrid scale," It In fact has some components In the 
resolved domain. A new modeling concept Involves exploiting this fact. 
We use the following decompositions : 


where 


ui + u| 


U^ = + u|' ' 


(C.1.1) 


u;(k) 


u^(k) ; 


u**»(k) i 




k > k. 


k^ Is the highest computed wavenumber. Therefore 

u"(k) = 11 u(k) 

G(k) 

We propose that 


Tfj « + C2m^j 


(C.l. 2) 


where 




(C.l. 3) 


and 


m 


Ij 



(C.l. 4) 
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? - 

Using the high wavenumber spectrum of Comte-Bellot and Corrsin (1971) 

D(k) » exp [^- I (C.1.5) 

We evaluate this expression for the highest computed wavenumber k^ 
using the experimental value, a = 1.3, to get 


E(k) = 3k exp 


r 3 

- j «(kn) 


E(k^) 


exp 


[- I 


where the Kolmogorov microscale n = (v/q3)lA, This gives Cj^ as 

■l/tl 


/: 


E(k) dk 


< T > 
ii 


(C.1.6) 


since m^^l^ = 0. To evaluate C2, we use Lumley's expression for the 
energy transfer spectrum : 


T(k) 


0.5tt ,5/2 

-372 

a 


E(0.4k) E^'^^(k) 


(C.1.7) 


and therefore 


T(k^) - < 

< "1 ”13.3 ^ 


This model is for homogeneous flow and therefore was not used in the 
present work. With more development, it may be used in flows with homo- 
geneity in two directions. .Jorge Bardina is testing models using simi- 
lar ideas and is getting very promising results. 
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Appendix D 
IMAGE FLOWS 


In solving problems in fluid mechanics , boundary conditions are 
often applied at finite boundaries. If this Is done for a flow which 
actually extends to Infinity, the boundary conditions imply artificial 
image flows. If one is not careful, the image flows may render the 
results meaningless. The new Fourier method described earlier puts the 
image flows infinitely far away, thus eliminating this difficulty. 

Suppose we want to compute the flow due to two vortices of the same 
sign. If we were to use discrete cosine and sine expansions on a finite 
uniform grid, we would be applying a no-stress boundary condition at 
some finite distance from the vortices. This implies an infinite array 
of pairs of vortices with alternating signs of vortlcity. In order to 
assess whether the image flows (a finite distance away) affect the com- 
puted solution, we performed the following numerical experiment. First 
we placed two pairs of vortices of opposite sign a distance d/2 above 
and below the x axis, as shown in Fig. 0.1. The vortices had ellipti- 
cal Gaussian distributions of vortlcity and a separation of distance 
c. The upper vortices will rotate about one another in a clockwise 
manner, while the lower pair will rotate counterclockwise. If the pairs 
are far enough apart not to affect each other, we should be able to 
shift the location of the lower pair by a distance c/2 in the x 
direction and get identical results. In performing this calculation, we 
used the mapping given by Eq. (3.3.5). We used the standard Fourier 
method for differentiating in the x direction and the second-order 
Adams-Bashforth scheme for the time advance. We used 16 grid points in 
the X direction and 128 in the z direction. The vortex pairs were 
centered at grid points 59 and 71 in z, and grid points 5 and 11 in 

X. The coefficient of the mapping function a was a = 192/tt (this 
gives Az = 1.5 near the origin); the x coordinates were Xj = 2(j- 
1); the dimensionless time step was selected so that (u^ ^^ /Ax + 
'^max^^^min^ At = .3. We did this calculation by solving the vortlcity 
equation : 
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3 2 

UW + Vtti = VV^U) 


(D.l) 


3 

Tt 


w + 




and the vector potential given by 


= - oj (D.2) 

to get 

3 I 3 I /•T^ Q\ 

u = -_t; w = (D.3) 

We computed until the dimensionless time T = = 1.0, at 

which point we compared the "turbulent" kinetic energy in the full do- 
main of the computation (with turbulence defined as the local deviation 
from an average in the x direction). For the case shown in Fig. D.l, 
the turbulent kinetic energy was the same at T = 1 and T = 0. How- 
ever, when the lower pair of vortices was shifted by c/2, the kinetic 
energy at T = 1 was double the energy at T = 0 . 

We thus conclude that a computation in a domain of height 1 .5 times 
the spacing between a pair of vortices would suffer tremendously from 
the influence of image flows . We also did the same calculation with 
d/c = 4 and found no significant image flow influence. 
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Fig. 1.3.1 The spatially developing mixing layer as 
created In the laboratory. 
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Fig. 4,1.1 Comparison of finite difference, new 
Fourier, and analytical solutions to the 
first-order convection equation. 
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Fig. 4.2.1 Comparison of the maximum dimensionless 
errors of finite-difference and new 
Fourier solutions to the heat or diffu- 
sion equation. 
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Fig. 5.3.1 


Mean velocity profile at various times 
for Case 6. 












Flg« 5.4.3. Normalized momentum thickness vs. Fig. 5.4.4. Dimensionless momentum thicknes growth 

dimensionless time for high“amplitude rate vs. normalized momentum thickness 

Cases 2 and 8 . for low-amplitude Cases 4 and 6 . 
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Fig. 5.4.5. Dimensionless momentum thickness growth 
rate vs. normalized momentum thickness 
for medixim- amplitude Cases 5, 7, 10, 11, 
and 9 . 



Fig. 5.4.6 Dimensionless momentum thickness growth 
rate vs . normalized momentxan thickness 
for high-amplitude Cases 1 and 8. 
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Fig. 5.4.7. Dimensionless momentum thickness growth 
rate vs . normalized momentum thickness 
for cases with seed #1. 
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Fig. 5.4.8. Dimensionless momentum thickness growth 
rate vs . normalized momentum thickness 
for cases with seed #2 . 
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Fig. 5.5.1. Dimensionless turbulent kinetic energy 
at the center of the layer vs. dimen- 
sionless time for low-amplitude Cases 4 
and 6. 


Fig. 5.5.2 R.M.S. values of the stteamwise velocity 
fluctuation at the center of the mixing 
layer versus ^ ~ ^ 

corresponds to T between 600 and 1200, 
depending on the initial boundary layer. 




Fig. 5.5.3 R.M.S. values of the gradient velocity 
fluctuation at the center of the mixing 
layer versus x/ r^ . See caption of Fig . 
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Fig. 5.5.5. Dimensionless turbulent kinetic energy 
at the center of the layer vs . normal- 
ized momentum thickness for medium- 
amplitude Cases 5, 7, 9, 10, and 11. 



Fig. 5.5.6. Dimensionless turbulent kinetic enefrgy 
at the center of the layer vs. normal- 
ized momentum thickness for high- 
amplitude Cases 2 and 8 . 
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5.7.1. 


Normalized Reynolds stress tensor pro- 
files for Case 6 at T = 0 . 


Fig. 5.7.2 


Normalized Reynolds stress tensor pro- 
files for Case 6 at T = 83. 

















Fig. 5.7-5. Normalized Reynolds stress tensor pro- 
files for Case 6 at T = 328. 


Fig . 5.7.6. 


Normalized Reynolds stress tensor pro 
files for Case 6 at T = 418. 







Fig. 5.7.7 


Normalized Reynolds stress tensor pro- 
files for Case 6 at T = 554. 


Fig. 5.7.8 


Normfitiized Reynolds stress tensor pro- 
files for Cases 5, 9, and 10 at T = 

0 . 0 . 







Fig. 5.7.9. Normalized Reynolds stress tensor pro- 
files for Case 7 at T = 0.0. 


Fig. 5.7.10. Normalized Reynolds stress tensor pro- 
files for Case 11 at T = 0.0. 












Fig. 5.7.13. Normalized Reynolds stress tensor pro- 
files for Case 9 at T = 66 . 


Fig. 5.7.14. Normalized Reynolds stress tensor prd* 
files for Case 11 at T = 67. 






Fig. 5.7.15. Normalized Reynolds stress tensor pro- Fig. 5.7.16. Normalized Reynolds stress tensor pro- 
files for Case 5 at T = 132. files for Case 7 at T = 132. 





Fig. 5.7.17. Normalized Reynolds stress tensqr pro' 
files for Case 9 at T = 131. 


Fig. 5.7.18. Normalized Reynolds stress tensor pro 
files for Case 11 at T = 131. 







Fig. 5.7.19. Normalized Reynolds stress tensor pro- 
files for Case 5 at T =» 197. 


Fig. 5.7.20. Normalized Reynolds stress tensor pro- 
files for Case 7 at T = 197. 










Fig. 5.7.23. Normalized Reynolds stress tensor pro- 
files for Case 5 at T = 295. 


Fig. 5.7.24. Normalized Reynolds stress tensor pro 
files for Case 7 at T = 295. 






Fig. 5.7.25. Normalized Reynolds stress tensor pro- 5.7.26. Normalized Reynolds stress tensor pro- 
files for Case 7 at T = 379. files for Case 2 at T = 33. 
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Fig. 5.8.1. Time history of the anisotropy tensor Fig. 5.8.2. Time history of the anisotropy tensor 

j for Case 4 . j for Case 6 . 






Fig. 5.8.3. Time history of the anisotropy tensor 

B . for Case 5 . 
ij 


Fig. 5.8.4. Time history of the anisotropy tensor 
j for Case 7 . 
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Fig. 5.8.5. Time history of the anisotropy tensor 
B. . for Case 9 . 


Fig. 5.8.6. Time history of the anisotropy tensor 

for Case 10. 






Fig. 5.8.7. Time history of the anisotropy tensor 

for Case 11. 


Fig. 5.8.8. Time history of the anisotropy tensor 

for Case 2 . 







Fig. 5.8.9. Time history of the anisotropy tensor 

for Case 8. 
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Fig. 5.9.1. Correlation coefficient at the center of 
the layer vs . normalized momenttna thick- 
ness for cases with seed #1 . 



Fig. 5.9.2. Correlation coefficient at the center of 
the layer vs . normalized momentum thick- 
ness for seed //2 . 
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Fig. 5.9.5 


Correlation coefficient at the center of 
the layer vs . normalized momentum thick- 
ness for high- amplitude Cases 2 and 8. 
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Fig. 5.10.1. Time history of alias-free accuracy and 
sample size parameters and rig. 
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Fig. 5.10.2. Time history of alias-free accuracy and 
sample size parameters t)^, ti^, and Hg. 






Fig. 5.10.3. Time history of alias-free accuracy and 


sample size parameters n 
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Fig. 5.10.4. Time history of alias-free accuracy and 
sample size parameters n^p, and Hg . 
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Fig. 5.10.5. Time history of alias-free accuracy and 
sample size parameters and Hg . 
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Fig. 5.11.2. Plan view of particle tracking grid for 
Case 6 at T = 161. 
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Fig. 5.11.3. Plan view of particle tracking grid for 
Case 6 at T = 200 . 


Fig. 5.11.4. Plan view of particle tracking grid for 
Case 6 at T = 280. 
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Fig. 5.11.5. Plan view of particle tracking grid for 5.11.6. Span view of particle tracking grid for 

Case 6 at T = 342 . Case 6 at T = 200 . 


Fig. 5.11.7. Span view of particle tracking grid for 
Case 6 at T = 342. 
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Fig. 5.11.9. Plan view of particle tracking grid for Fig. 5.11.10. Plan view of particle tracking grid for 
Case 8 at T = 76. Case 10 at T = 76. 
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Fig. 5.11.11. Plan view of particle tracking grid for 
Case 10 at T = 97. 


Fig. 5.11.12. Plan view of particle tracking grid for 
Case 10 at T = 151. 
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Fig. 5.11.13. Plan view of particle tracking grid for Fig. 5.11.14. Plan view of particle tracking grid for 
Case 10 at T = 171. Case 10 at T = 231. 





Fig. 5.11.17. Plan view of particle tracking grid for 
Case 9 at T = 220. 


Fig. 5.11.18. Span view of particle tracking grid for 
Case 9 at T = 220 . 
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Fig. 5.11.19. Plan view of particle tracking grid for Fig. 5.11.20. Plan view of particle tracking grid for 

Case 11 at T = 97. Case 11 at T = 151. 


123 


Fig. 5.11.21. Plan view of particle tracking grid for 
Case 11 at T = 231. 




Fig. 5.11.22 Span view of particle tracking grid for 
Case 11 at T = 171. 
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• 5.11.23. Span view of particle tracking grid for 
Case 11 at T = 231. 




Fig. D.l Schematic view of vortex pairing, Image 
flow study . 


125 


End of Document 



